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)
becomesfoo_(&i,&x)
for a functionextern "C" foo_(int*, float*)
that takes pointer arguments. C++ allows for the alternative declarationextern "C" foo_(int&, float&);
, which can be invoked with the more natrual syntaxfoo_(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’sA[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>
orcomplex<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/Accelerate.framework/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\times 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 int integer;
// or typedef long int integer; on some systems
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.