How to Write an HTML5 App? A tutorial about writing HTML5 Apps, using HTML, JavaScript and CSS. … read more

Calculating the Square Root of a Matrix. How to use LAPACK in a C program to calculate the square root of a positive definite matrix. … read more

Matrix Analysis and Algorithms. Lecture Notes: Numerical Linear Algebra … read more

By Jochen Voss, last updated 2012-02-18

This page explains how to use numerical linear algebra packages of the BLAS/LAPACK family in C programs on a Debian GNU/Linux system. My goal is to efficiently solve systems of linear equations. Any feedback about this page is very welcome.

- Overview
- Installation
- Using BLAS
- Using LAPACK the Easy Way
- Using LAPACK the Hard Way
- Outlook
- References

There is a jungle of different packages and standards and there are implementations in both C and Fortran. The more important packages are presented in the following list.

- BLAS
- The
*Basic Linear Algebra Subprograms*are high qualitybuilding block

routines for performing basic vector and matrix operations. Level 1 BLAS do vector-vector operations, Level 2 BLAS do matrix-vector operations, and Level 3 BLAS do matrix-matrix operations. There seem to be many different implementations of these standardised routines.Originally the BLAS library is a Fortran library. As we will see below, it is possible to call the functions from this library from a C program. There is also a standardised C language interface, named

`cblas`

, for the library. - LINPACK
- LINPACK is a collection of Fortran subroutines that analyse and solve linear equations and linear least-squares problems. LINPACK seems to be superseded by LAPACK.
- EISPACK
- EISPACK is a collection of Fortran subroutines that compute the eigenvalues and eigenvectors of matrices. EISPACK seems to be superseded by LAPACK.
- LAPACK
- LAPACK provides routines for solving systems of simultaneous linear
equations, least-squares solutions of linear systems of equations,
eigenvalue problems, and singular value problems. The associated matrix
factorisations (LU, Cholesky, QR, SVD, Schur, generalised Schur) are also
provided, as are related computations such as reordering of the Schur
factorisations and estimating condition numbers. Dense and banded matrices
are handled, but not general sparse matrices. In all areas, similar
functionality is provided for real and complex matrices, in both single
and double precision. This is the package we are going to use.
Originally the LAPACK library is a Fortran library but it is possible to call the functions from a C program. There is also a non-standard, partial C language interface for this libary, named

`clapack`

. - ATLAS
- The
*Automatically Tuned Linear Algebra System*is a fast implementation of the BLAS and of a subset of LAPACK. This implementation of LAPACK is commonly used on Linux.

On my Debian GNU/Linux system I installed the following packages to make use of the linear algebra routines. On non-Debian systems the package names might be different.

`atlas3-sse2-dev`

- This is the ATLAS version of the BLAS and LAPACK libraries. Because my machine has a Pentium IV processor I chose the SSE2 variant of the package.
`lapack3-dev`

- This contains the libatlas.so and (per dependency) libatlas.so symlinks.
`refblas3-dev`

- This contains the BLAS header file cblas.h.
`refblas3-doc`

- This documents the BLAS routines. The relevant files are blas2-paper.ps.gz, blas3-paper.ps.gz. These describe the semantics and the Fortran interface. You also need to read cblas.ps.gz (see below) to find the corresponding C function names.
`lapack3-doc`

- This package contains a copy of the LAPACK user's guide. Point your browser to file:///usr/share/doc/lapack3-doc/lug/index.html to access it. Also included is the LAPACK quick reference card lapackqref.ps.gz.
`atlas3-doc`

- This contains mostly documentation about the automatic tuning process, which is irrelevant for our purposes. But in /usr/share/doc/atlas3-doc/ it also contains a quick reference sheet cblasqref.ps.gz and the reference document cblas.ps.gz for the BLAS C bindings.

To test the BLAS routines we want to perform a simple matrix-vector
multiplication. Reading the file blas2-paper.ps.gz
we find that the name of the corresponding Fortran function is
`DGEMV`

. The text blas2-paper.ps.gz
also explains the meaning of the arguments to this function. In
cblas.ps.gz we find that the corresponding C
function name is `cblas_dgemv`

. The following example uses this
function to calculate the matrix-vector product

/ 3 1 3 \ / -1 \
| 1 5 9 | * | -1 |.
\ 2 6 5 / \ 1 /

Example file testblas.c:

#include <stdio.h> #include <cblas.h> double m[] = { 3, 1, 3, 1, 5, 9, 2, 6, 5 }; double x[] = { -1, -1, 1 }; double y[] = { 0, 0, 0 }; int main() { int i, j; for (i=0; i<3; ++i) { for (j=0; j<3; ++j) printf("%5.1f", m[i*3+j]); putchar('\n'); } cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1.0, m, 3, x, 1, 0.0, y, 1); for (i=0; i<3; ++i) printf("%5.1f\n", y[i]); return 0; }

To compile this program we use the following command.

cc testblas.c -o testblas -lblas -lm

The output of this test program is

3.0 1.0 3.0 1.0 5.0 9.0 2.0 6.0 5.0 -1.0 3.0 -3.0

which shows that everything worked fine and that we did not even use the transposed matrix by mistake.

To test the LAPACK routines we will solve a simple system of linear
equations. In the
linear equations
section of the LAPACK user's guide we find that the corresponding Fortran
function name is `DGESV`

.

The Atlas library provides a C wrapper for the Fortran function
`DGESV`

which is named `clapack_dgesv`

. You can find
out about the available C wrapper functions by reading the
lapackqref.ps.gz quick reference page, which is
found in the /usr/share/doc/atlas3-doc/ directory.

To compile the following example, you need the `clapack.h`

header file from the Atlas distribution. Since revision
`3.6.0-6`

this file can be found in the
`atlas3-headers`

Debian package.

Example file testlapack.c:

#include <stdio.h> #include <atlas_enum.h> #include "clapack.h" double m[] = { 3, 1, 3, 1, 5, 9, 2, 6, 5 }; double x[] = { -1, 3, -3 }; int main() { int ipiv[3]; int i, j; int info; for (i=0; i<3; ++i) { for (j=0; j<3; ++j) printf("%5.1f", m[i*3+j]); putchar('\n'); } info = clapack_dgesv(CblasRowMajor, 3, 1, m, 3, ipiv, x, 3); if (info != 0) fprintf(stderr, "failure with error %d\n", info); for (i=0; i<3; ++i) printf("%5.1f %3d\n", x[i], ipiv[i]); return 0; }

To compile this program we use the following command.

cc testlapack.c -o testlapack -llapack_atlas -llapack -lblas -latlas -lm

The resulting program obviously depends on the Atlas libraries and won't run with any other BLAS implementation.

The program reverses the matrix-vector multiplication from our first
example and correctly recovers the vector `c`

. The output of
this test program is

3.0 1.0 3.0 1.0 5.0 9.0 2.0 6.0 5.0 -1.0 0 -1.0 2 1.0 2

As a more interesting test for the LAPACK routines we solve a
tridiagonal system of linear equations. The Fortran function name as found
in the
linear equations
section of the LAPACK user's guide is `DGTSV`

.

The lapackqref.ps.gz quick reference page reveals that Atlas does not provide a C wrapper for this function. Thus we have to improvise one on our own.

The following example uses this method to solve the system

/ 2 -1 0 0 0 \ / x0 \ / 1 \
| -1 2 -1 0 0 | | x1 | | 2 |
| 0 -2 3 -1 0 | * | x2 | = | 3 |.
| 0 0 -1 3 -2 | | x3 | | 2 |
\ 0 0 0 -1 1 / \ x4 / \ 1 /

From the
Tridiagonal and Bidiagonal Matrices
page in the LAPACK documentation we know that this tridiagonal matrix
should be stored in three arrays of length 4 (sub-diagonal), 5 (diagonal),
and 4 (super-diagonal). These are one-dimensional arrays so we won't get
confused by the fact that Fortran expects two-dimensional arrays to be in
column-major order. We can access the function `DGTSV`

from the
`lapack`

library under the name `dgtsv_`

, the full
argument list can be deduced either from the initial comment of the
dgtsv.f Fortran
file in the netlib archive or from the
lapackqref.ps.gz quick reference card.

Example file testdgtsv.c:

#include <stdio.h> double l[] = { -1, -2, -1, -1 }; double d[] = { 2, 2, 3, 3, 1 }; double u[] = { -1, -1, -1, -2 }; double x[] = { 1, 2, 3, 2, 1 }; static long dgtsv(long N, long NRHS, double *DL, double *D, double *DU, double *B, long LDB) { extern void dgtsv_(const long *Np, const long *NRHSp, double *DL, double *D, double *DU, double *B, const long *LDBp, long *INFOp); long info; dgtsv_(&N, &NRHS, DL, D, DU, B, &LDB, &info); return info; } int main() { int i, info; info = dgtsv(5, 1, l, d, u, x, 5); if (info != 0) fprintf(stderr, "failure with error %d\n", info); for (i=0; i<5; ++i) printf("%5.1f\n", x[i]); return 0; }

To compile this program we use the following command.

cc testdgtsv.c -o testdgtsv -llapack -lblas

This time the program is completely independent of the Atlas libraries. If the optimised Atlas versions are present it will use them, if they are absent it will just use the reference implementation instead.

The output of the test program is

6.5 12.0 15.5 19.5 20.5

which is indeed the solution to the above system of equations.

Remark. The same source code can be compiled under MacOS X, but the linker command has to be changed to

cc testdgtsv.c -o testdgtsv -framework veclib

Until now we have illustrated the basics of how to make use of BLAS and LAPACK library functions from within C programs. A more exhaustive example, illustrating the Fortran conventions for storing two-dimensional arrays in memory, is presented on my page about calculating the square root of a matrix.

- the ATLAS homepage
- the LAPACK homepage, containing the LAPACK user's guide
- the BLAS homepage
- my numerical linear algebra lecture notes.

Copyright © 2012, Jochen Voss. All content on this website (including text, pictures, and any other original works), unless otherwise noted, is licensed under a Creative Commons Attribution-Share Alike 3.0 License.