Nicolas LIMARE / pro / notes / 2014 / 10 / Using BLAS and LAPACK from C/C++

How can we call the BLAS and LAPACK libraries from a C code without being tied to an implementation? For BLAS, there is CBLAS, a native C interface. For LAPACK, the native C interface is LAPACKE, not CLAPACK. If you don't have LAPACKE, use extern Fortran declarations.

BLAS and LAPACK

The reference Fortran code for BLAS and LAPACK defines de facto a Fortran API, implemented by multiple vendors with code tuned to get the best performance on a given hardware.

In fact, since libraries are often distributed in binary form (think Linux distro packages for example), BLAS and LAPACK could be written in any language as long as that a Fortran code calling a BLAS or LAPACK subroutine with the syntax used on Netlib can compile and link with the library and the resulting executable produces the same result the Netlib code.

Well, I have not touched Fortran for years bus as long as I remember Fortran binary interface (ABI) is not standardized. An object file produced by a compiler may not be usable by another compiler (cf g95 vs. gfortran vs. ifort). But let's just suppose here that we are only working with one compiler with a well defined ABI. So, if we know the binary calling convention for a Fortran function, the the Fortran API also defines an ABI.

On the C side, the C ABI is standardized: from a source-level function definition, we know how the binary-level call will happen. Moreover, there is a C type matching every Fortran scalar type used in BLAS and LAPACK. As a result, the Fortran source API for BLAS/LAPACK plus assumptions about the Fortran compiler result in a C source API for BLAS/LAPACK.

For example, we can access the BLAS dgemm and LAPACK dgesv subroutines ...

SUBROUTINE dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DOUBLE PRECISION alpha,beta
INTEGER k,lda,ldb,ldc,m,n
CHARACTER transa,transb
DOUBLE PRECISION a(lda,*),b(ldb,*),c(ldc,*)

SUBROUTINE dgesv(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
INTEGER info, lda, ldb, n, nrhs
INTEGER ipiv(*)
DOUBLE PRECISION a(lda, *), b(ldb, *)

... as the external C functions dgemm_ and dgesv_ ...

extern dgemm_(char * transa, char * transb, int * m, int * n, int * k,
              double * alpha, double * A, int * lda,
              double * B, int * ldb, double * beta,
              double *, int * ldc);

extern dgesv_(int * n, int * nrhs, double * A, int * lda,
              int * ipiv, double * B, int * ldb, int * info);

... and this stands for any BLAS/LAPACK subroutine.

Note that since this Fortran API to C ABI is compiler-dependent, no BLAS library provides a lapack.h header with all these extern declarations. You have to do it yourself, and trust the compilers for long term validity. Or you may prefer a truly portable, standard interface.

CBLAS

This interface exists for BLAS, it's called CBLAS. Following on the dgemm example, we now have this new C API/ABI:

void cblas_dgemm(const enum CBLAS_ORDER Order,
                 const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_TRANSPOSE TransB,
                 const int M, const int N, const int K,
                 const double alpha, const double  *A, const int lda,
                 const double *B, const int ldb, const double beta,
                 double *C, const int ldc)

CBLAS is in fact just a very thin layer (no performance penalty) over BLAS, with these benefits:

  • standard and simple, implemented by every BLAS library;
  • compiler-independent;
  • scalars passed by value, arrays by pointers;
  • const to protect read-only parameters;
  • enum and macros instead of characters with "magic values";
  • Order parameter to use row-major or column-major arrays without juggling with transpositions.

If you use BLAS from a C code, there is no reason to use the Fortran interface directly. Just #include <cblas.h> (provided with every BLAS library) and link with your preferred BLAS.

CLAPACK?

There also is something on Netlib called CLAPACK, but it is not a standard C interface to the Fortran LAPACK. CLAPACK is a translation of the Fortran LAPACK code into C.

Still on the dgesv example, now in CLAPACK it is defined as ...

int dgesv_(integer *n, integer *nrhs, doublereal *a,
           integer *lda, integer *ipiv, doublereal *b,
           integer *ldb, integer *info);

... with integer and doublereal typedef'd as int and double somewhere else. Based on 100% C code, this interface is compiler-independent. But that's all you can get: it is not implemented by other vendors and does not provide any of the other benefits of a native C interface like CBLAS does.

Finally, CLAPACK is mentioned as "no longer maintained" on Netlib, so let's just avoid this trap.

There is also something called clapack.h in ATLAS, with functions like ...

int clapack_dgesv(const enum CBLAS_ORDER Order,
                  const int N, const int NRHS,
                  void *A, const int lda, int *ipiv,
                  void *B, const int ldb);

This something else, a C interface provided by ATLAS for some of the LAPACK routines in ATLAS, entirely different from Netlib CLAPACK. And ATLAS doesn't provide every LAPACK subroutine. So it's time to conclude on the absence of anything standard and portable called CLAPACK. As written by Clint Whaley (author of ATLAS), "there has been no official standardization of a C interface to LAPACK ... your best bet is to call the Fortran77 interface".

Well that was sadly true until ...

LAPACKE

Here it is. LAPACKE is relatively recent (proposed in 2008, accepted in 2010) but it's the closest we have to a native C source-level interface to the original Fortran API. LAPACKE is to LAPACK what CBLAS is to BLAS, with this interface:

lapack_int LAPACKE_dgesv(int matrix_order, lapack_int n, lapack_int nrhs,
                         double * a, lapack_int lda, lapack_int * ipiv,
                         double * b, lapack_int ldb);

This one is standardized, as mentioned on Netlib: "Standard C language APIs for LAPACK". Moreover, as LAPACKE is a thin layer over LAPACK, which itself calls BLAS to perform the real expensive computations, I think that it is possible to use the liblapacke.so library from Netlib with any LAPACK and BLAS back-end. Everything you need to the good linking options, in the correct order. Finally, LAPACKE is also included in the MKL package.

All Together

The BLAS case is clear: we must use the CBLAS interface, available with every BLAS library. For LAPACK, we should use LAPACKE if we have it. If not, we could either use the LAPACKE layer from Netlib, and possibly the LAPACK layer too, over our preferred BLAS library, or directly link to the Fortran interface if we know the compiler ABI.

For reference, here is the list of Fortran and C interfaces available (ignoring unmaintained CLAPACK) in major BLAS/LAPACK implementations, with the package, shared library and C header names in Debian 8.0 Jessie:

Netlib (pkg libblas-dev, liblapack-dev, liblapacke-dev)

  • BLAS (/usr/lib/libblas/libblas.so),
  • CBLAS (/usr/lib/libblas/libblas.so, /usr/include/cblas.h)
  • LAPACK (/usr/lib/lapack/liblapack.so)
  • LAPACKE (/usr/lib/liblapacke.so, /usr/include/lapacke.h)

Atlas (pkg libatlas-dev)

  • BLAS (/usr/lib/atlas-base/libf77blas.so)
  • CBLAS (/usr/lib/atlas-base/libcblas.so, /usr/include/atlas/cblas.h)
  • LAPACK (/usr/lib/atlas-base/liblapack_atlas.so, incomplete)

OpenBLAS (pkg libopenblas-dev)

  • BLAS (/usr/lib/openblas-base/libblas.so)
  • CBLAS (/usr/lib/openblas-base/libblas.so, /usr/include/openblas/cblas.h)
  • LAPACK (/usr/lib/openblas-base/liblapack.so)

GSL (pkg libgsl0-dev)

  • CBLAS (/usr/lib/libgslcblas.so, /usr/include/gsl/gsl_cblas.h)

MKL

  • BLAS ($MKLROOT/lib/intel64/libmkl_<C>_<P>.so)
  • CBLAS, ($MKLROOT/lib/intel64/libmkl_<C>_<P>.so, $MKLROOT/include/mkl_cblas.h)
  • LAPACK ($MKLROOT/lib/intel64/libmkl_<C>_<P>.so)
  • LAPACKE ($MKLROOT/lib/intel64/libmkl_<C>_<P>.so, $MKLROOT/include/mkl_lapacke.h)

<C> is the compiler ABI (gf for GFortran and intel for IFort) and <P> is the 64-bit model, lp64 or ilp64.