Numerical Linear Algebra Packages on Linux

By Jochen Voss, last updated February 18, 2012.

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

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 quality building 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.

Installation

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.

Using BLAS

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.

Using LAPACK the Easy Way

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

Using LAPACK the Hard Way

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

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

Outlook

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.

References

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

Numerical Linear Algebra Packages on Linux

Contents

Overview

Installation

Using BLAS

Using LAPACK the Easy Way

Using LAPACK the Hard Way

Outlook

References