Understanding Lapack Calls in C++ with a Simple Example

Understanding LAPACK calls in C++ with a simple example

You need to factor the matrix (by calling dgetrf) before you can solve the system using dgetrs. Alternatively, you can use the dgesv routine, which does both steps for you.

By the way, you don't need to declare the interfaces yourself, they are in the Accelerate headers:

// LAPACK test code

#include <iostream>
#include <vector>
#include <Accelerate/Accelerate.h>

using namespace std;

int main()
{
    char trans = 'N';
    int dim = 2;    
    int nrhs = 1;
    int LDA = dim;
    int LDB = dim;
    int info;

    vector<double> a, b;

    a.push_back(1);
    a.push_back(1);
    a.push_back(1);
    a.push_back(-1);

    b.push_back(2);
    b.push_back(0);

    int ipiv[3];

    dgetrf_(&dim, &dim, &*a.begin(), &LDA, ipiv, &info);
    dgetrs_(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);

    std::cout << "solution is:";    
    std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
    std::cout << "Info = " << info << std::endl; 

    return(0);
}

call lapack and blas in c++

Vendor supplied BLAS and LAPACK implementations usually include symbols with the trailing underscore, because that is how Fortran 77 compilers worked originally. The modern gfortran behavior is to also add the trailing underscore for compatibility, but there is the -fno-underscoring option to turn it off.

For code you compile yourself with gfortran such as the reference BLAS and LAPACK, you can choose whether to return complex values directly or by using the indirect result pointer argument. To get the indirect behavior, compile with -ff2c. The default gfortran behavior is to return complex values directly.

The simplest way to avoid the macros and wraps in your code is to assume trailing underscores in the names and indirect return values for complex results in an added first argument. That will be compatible with vendor libraries. Then compile BLAS and LAPACK with -ff2c to make it work there.

For maximum flexibility you can use wrap functions. Within the wrapper you should only need to worry about whether complex arguments are returned directly or not, you should not need special handling for each different library.

complex<double> zdotc_wrap(int n, const complex<double>*x, int incx, const complex<double>*y, int incy) 
{
  #if defined(FORTRAN_COMPLEX_FUNCTIONS_RETURN_VOID)
    complex<double> result;
    zdotc_(&result, &n, x, &incx, y, &incy);
    return result;
  #else
    return zdotc_(&n, x, &incx, y, &incy);
  #endif
}

In BLAS there are only a handful of functions that need to be wrapped: CDOTU CDOTC ZDOTU ZDOTC. In LAPACK there are just CLADIV ZLADIV (I think).

Link to Fortran library (Lapack) from C++

As you have already discovered, hacks like adding the underscore are platform-specific by their nature. Supporting a big range of platforms and compilers that way "by hand" requires a large amount of unpleasant and tedious work.

The easiest way to get portability is to use LAPACKE, the official C interface to LAPACK. This also has the added benefit of saving you from having to re-declare all the functions you need.

Here's a simple example:

#include <iostream>
#include <lapacke.h>    

int main()
{
    // By using lapack_int, we also support LAPACK-ILP64
    lapack_int major = 0;
    lapack_int minor = 0;
    lapack_int patch = 0;
    LAPACKE_ilaver(&major, &minor, &patch);
    std::cout << major << "." << minor << "." << patch << std::endl;
    return 0;
}

More information can be found in the official documentation.

Note that LAPACKE only handles LAPACK, if you also need BLAS routines, you can get those from CBLAS.

Computing the inverse of a matrix using lapack in C

First, M has to be a two-dimensional array, like double M[3][3]. Your array is, mathematically speaking, a 1x9 vector, which is not invertible.

• N is a pointer to an int for the
  order of the matrix - in this case,
  N=3.

• A is a pointer to the LU
  factorization of the matrix, which
  you can get by running the LAPACK
  routine dgetrf.

• LDA is an integer for the "leading
  element" of the matrix, which lets
  you pick out a subset of a bigger
  matrix if you want to just invert a
  little piece. If you want to invert
  the whole matrix, LDA should just be
  equal to N.

• IPIV is the pivot indices of the
  matrix, in other words, it's a list
  of instructions of what rows to swap
  in order to invert the matrix. IPIV
  should be generated by the LAPACK
  routine dgetrf.

• LWORK and WORK are the "workspaces"
  used by LAPACK. If you are inverting
  the whole matrix, LWORK should be an
  int equal to N^2, and WORK should be
  a double array with LWORK elements.

• INFO is just a status variable to
  tell you whether the operation
  completed successfully. Since not all
  matrices are invertible, I would
  recommend that you send this to some
  sort of error-checking system. INFO=0 for successful operation, INFO=-i if the i'th argument had an incorrect input value, and INFO > 0 if the matrix is not invertible.

So, for your code, I would do something like this:

int main(){

    double M[3][3] = { {1 , 2 , 3},
                       {4 , 5 , 6},
                       {7 , 8 , 9}}
    double pivotArray[3]; //since our matrix has three rows
    int errorHandler;
    double lapackWorkspace[9];

    // dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix 
    // called A, sending the pivot indices to IPIV, and spitting error 
    // information to INFO.
    // also don't forget (like I did) that when you pass a two-dimensional array
    // to a function you need to specify the number of "rows"
    dgetrf_(3,3,M[3][],3,pivotArray[3],&errorHandler);
    //some sort of error check

    dgetri_(3,M[3][],3,pivotArray[3],9,lapackWorkspace,&errorHandler);
    //another error check

    }

Matrix operations in C++ (using Blas/Lapack or some other alternative)

I assume that if you are new to C++, you are also new to C and Fortran. In that case I would definitely suggest to you, not to start with BLAS/LAPACK, at least not without a nice C++ wrapper.

My suggestion would be to have a look at Eigen which offers a much easier start to matrix operations using native C++ code. You can have a look at their tutorial to get started. The Eigen performance is said to be comparable to that of BLAS/LAPACK. See e.g. their benchmark. However I didn't test that myself.

If you really want to go low level and use BLAS/LAPACK, have a look at the available functions of cBlas (the C-Wrapper of BLAS) and LAPACK. Additionally, you can find some examples how to use Lapacke (The C-Wrapper of LAPACK) here. But don't expect things to be nice and well documented!

To finally give an answer to your question: Here is a code snipped I wrote some time ago for benchmarking. The code creates two random matrices A and B and multiplies them into the matrix C.

#include <random>
#include <cblas.h>

int main ( int argc, char* argv[] ) {

    // Random numbers
    std::mt19937_64 rnd;
    std::uniform_real_distribution<double> doubleDist(0, 1);

    // Create arrays that represent the matrices A,B,C
    const int n = 20;
    double*  A = new double[n*n];
    double*  B = new double[n*n];
    double*  C = new double[n*n];

    // Fill A and B with random numbers
    for(uint i =0; i <n; i++){
        for(uint j=0; j<n; j++){
            A[i*n+j] = doubleDist(rnd);
            B[i*n+j] = doubleDist(rnd);
        }
    }

    // Calculate A*B=C
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0, A, n, B, n, 0.0, C, n);

    // Clean up
    delete[] A;
    delete[] B;
    delete[] C;

    return 0;
}

---

Related Topics