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.
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.
building blockroutines 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.
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.
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
lapack3-dev
refblas3-dev
refblas3-doc
lapack3-doc
atlas3-doc
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
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
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.
Copyright © 2012 Jochen Voss. All content on this website (including text, pictures, and any other original works), unless otherwise noted, is licensed under the CC BY-SA 4.0 license.