Linking to the LAPACK libraries from C++

What is LAPACK?

LAPACK is a collection of high-performance linear algebra routines written in FORTRAN and built on top of BLAS. LAPACK is usually provided as a shared library compiled either from the FORTRAN source or from C code automatically generated by the f2c conversion utility. The latter form is referred to as CLAPACK.

Calling conventions

In either case, the libraries obey FORTRAN calling conventions, so there are several things to keep in mind:

• FORTRAN subroutines correspond to C functions that return void.
• In C++ code, the function prototypes should be declared extern "C" in order to turn off name mangling.
• Although FORTRAN is case insensitive, most compilers adopt the convention of making routines available to the linker in all lower-case letters. The name is sometimes followed by an underscore. For example, using g77 and gcc, the FORTRAN subroutine DSPEVD is accessed from C/C++ as dspevd_. Since the presence of the trailing underscore is compiler-dependent, it may be best to define

  #ifdef FORTRAN_TRAILING_UNDERSCORE
  #define F77NAME(x) x##_
  #else
  #define F77NAME(x) x
  #endif
  

and refer to F77NAME(dspevd) throughout your code.
• FORTRAN uses pass by reference semantics, so call FOO(i,x) becomes foo_(&i,&x) for a function extern "C" foo_(int*, float*) that takes pointer arguments. C++ allows for the alternative declaration extern "C" foo_(int&, float&);, which can be invoked with the more natrual syntax foo_(i,x).
• A key difference for the purposes of linear algebra calculations is that C/C++ store matrices (two-dimensional C arrays) in row-major order indexed from zero, whereas FORTRAN stores them in column-major order indexed from one. In other words, FORTRAN's A(3,7) is C's A[6][2].
• There is not always an exact correspondence between FORTRAN and C/C++ data types. The main culprits are the INTEGER and LOGICAL types, which are ints and unsigned ints that may or may not correspond to a long, and the COMPLEX and DOUBLE COMPLEX types, which are C structures and must be explicitly cast to a complex<float> or complex<double>. Note that CLAPACK distributions are bundled with a header file clapack.h that defines function prototypes and typedefs corresponding FORTRAN's numerical types. On Mac OS X, there is a file /System/Library/Frameworks/vecLib.framework/Headers/clapack.h:

  #ifndef __CLAPACK_H
  #define __CLAPACK_H
  
  #ifdef __cplusplus
  extern "C" {
  #endif
  	
  #if defined(__LP64__) /* In LP64 match sizes with the 32 bit ABI */
  typedef int 		__CLPK_integer;
  typedef int 		__CLPK_logical;
  typedef float 		__CLPK_real;
  typedef double 		__CLPK_doublereal;
  typedef __CLPK_logical 	(*__CLPK_L_fp)();
  typedef int 		__CLPK_ftnlen;
  #else
  typedef long int 	__CLPK_integer;
  typedef long int 	__CLPK_logical;
  typedef float 		__CLPK_real;
  typedef double 		__CLPK_doublereal;
  typedef __CLPK_logical 	(*__CLPK_L_fp)();
  typedef long int 	__CLPK_ftnlen;
  #endif
  	
  typedef struct { __CLPK_real r, i; } __CLPK_complex;
  typedef struct { __CLPK_doublereal r, i; } __CLPK_doublecomplex;
  
  // + function prototypes
  
  #ifdef __cplusplus
  }
  #endif
  #endif /* __CLAPACK_H */
  

Example code

Here's a program that diagonalizes a real, symmetric N×N matrix arranged in the upper-triangular packed storage format.

#ifdef __APPLE__
#include <Accelerate/Accelerate.h>
#include <vecLib/clapack.h>
typedef __CLPK_doublereal doublereal;
typedef __CLPK_integer integer;
#else
typedef double doublereal;
typedef long int integer;
extern "C" void F77NAME(dspevd)
   (char*, char*, integer*, doublereal*, doublereal*,
    doublereal*, integer*, doublereal*, integer*,
    integer*, integer*, integer*);
#endif

integer eigensolve(vector<doublereal> &H, 
                   vector<doublereal> &Eval, 
                   vector<doublereal> &Evec)
{
   // Solve the eigenvalue problem with LAPACK's dsepvd routine

   integer N = Eval.size();
   assert(H.size() == N*(N+1)/2);
   assert(Evec.size() == N*N);

   integer info;
   char jobz='V';
   char uplo='U';
   vector<doublereal> work(1+6*N+N*N);
   integer lwork = work.size();
   vector<integer> iwork(3+5*N);
   integer liwork = iwork.size();

   F77NAME(dspevd)(&jobz,&uplo,&N,&(H[0]),&(Eval[0]),&(Evec[0]),&N,
           &(work[0]),&lwork,&(iwork[0]),&liwork,&info);

   return info;
}

Compiling and linking

On Mac OS X systems, the LAPACK libraries are provided as a part of the Accelerate framework. Linux requires that you link to /usr/lib/liblapack.so, /usr/lib/libblas.so, and /usr/lib/libg2c.so.

$ more makefile
CC=g++
CFLAGS=-O2 -ansi -pedantic

## for Mac OS X
LIBS=-framework Accelerate

## For Linux
#LIBS=-llapack -lblas -lg2c -lm

example: example.cpp
        $(CC) -o $@ $@.cpp $(CFLAGS) $(LIBS)

clean:
        rm example 2>&-