Feeds:
Posts
Comments

Archive for the ‘lapack’ Category

to call a LAPACK routine from your C++ code i used to do the following

extern "C"
{
  //http://www.netlib.org/lapack/double/dgesv.f
  void dgesv_(int* n, int *nrhs, double* a, int* lda, int* ipiv, double* b, int* ldb, int* info);
}

int main()
{
  int n = 2;
  int nrhs = 2;
  double a[] = {1,3,2,4}; //colum-major mode 
  int lda = n;
  int* ipiv = new int[n];
  double b[] = {19, 43, 22, 50}; //column-major mode
  int ldb = n;
  int info;
  
  dgesv_(&n, &nrhs, a, &lda, ipiv, b, &ldb, &info);

  delete [] ipiv;
  return 0;
}

the above code solves the system of linear equations \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}\mathbf{X} = \begin{bmatrix} 19 & 22 \\ 43 & 50 \end{bmatrix}. notice that in the code i had to store the matrices in column-major mode that is used by FORTRAN. i generally write helper routines to convert from row-major mode to column major mode before the subroutine call and then convert the results back to row-major mode after returning from the call. also one has to wrap the call to the subroutine in an extern “C” wrapper. both of these were annoyances that more or less kept me away from LAPACK, i could get most of my work done using GSL. while GSL provides a BLAS interface it does not provide any LAPACK interface and i have always wished to be able to use LAPACK.

it turns out that the guys at the intel math kernel library(imkl) and lapack guys collaborated to bring a C language interface to lapack . on archlinux i obtained the lapacke package from the AUR. below is the code using lapacke

#include <lapacke.h> // <--- header file
int main()
{
  double a[] = {1,2,3,4}; //NO need for column-major mode
  double b[] = {19, 22, 43, 50}; //NO need for column-major mode
  int n = 2;
  int nrhs = 2;
  int lda = n;
  int ipiv[n];
  int ldb = n;

  int info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, a, lda, ipiv, b, ldb);
 
  return 0;
}

compile the above code using g++ main.cpp -llapacke. notice that by using LAPACK_ROW_MAJOR on LAPACK_COL_MAJOR one can specify the storage mode for matrices and therefore there is no need to preprocess the matrices before calling the routine.

i am aware of the CLAPACK project as well which is an f2c’d version of the fortran version but i always wanted a pure C version. this i thought would fill my needs. alas that is NOT the case. i will not be using lapacke because it is not c++0x compatible 😦 and i use a lot of c++0x. i wish GSL would provide a lapack interface as well! my search for a scientific computing library for C++ continues.

Read Full Post »