Understanding LAPACK calls in C++ with a simple example

Question

I am a beginner with LAPACK and C++/Fortran interfacing. I need to solve linear equations and eigenvalues problems using LAPACK/BLAS on Mac OS-X Lion. OS-X Lion provides optimized BLAS and LAPACK libraries (in /usr/lib) and I am linking these libraries instead of downloading them from netlib.

My program (reproduced below) is compiling and running fine, but it is giving me wrong answers. I have researched in the web and Stackoverflow and the issue may have to deal with how C++ and Fortran store arrays in differing formats (row major vs Column major). However, as you will see in my example, the simple array for my example should look identical in C++ and fortran. Anyway here goes.

Lets say we want to solve the following linear system:

x + y = 2

x - y = 0

The solution is (x,y) = (1,1). Now I tried to solve this using Lapack as follows

// LAPACK test code

#include<iostream>
#include<vector>

using namespace std;
extern "C" void dgetrs(char *TRANS, int *N, int *NRHS, double *A, 
                      int *LDA, int *IPIV, double *B, int *LDB, int *INFO );

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];


    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);
}

This code was compiled as follows:

g++ -Wall -llapack -lblas lapacktest.cpp

On running this, however, I get the solution as (-2,2) which is obviously wrong. I have tried all combination of row/column re-arrangement of my matrix a. Also observe the matrix representation of a should be identical in row and column formats. I think getting this simple example to work will get me started with LAPACK and any help will be appreciated.

Solution 1

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);
}

Solution 2

For those who don't want bother with the Accelerate Framework, I provide the code of Stephen Canon (thanks to him, of course) with nothing but pure library linking

// LAPACK test code
//compile with: g++ main.cpp -llapack -lblas -o testprog

#include <iostream>
#include <vector>

using namespace std;

extern "C" void dgetrf_(int* dim1, int* dim2, double* a, int* lda, int* ipiv, int* info);
extern "C" void dgetrs_(char *TRANS, int *N, int *NRHS, double *A, int *LDA, int *IPIV, double *B, int *LDB, int *INFO );

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);
}

And about the manual, there's a full PDF version available at Intel's website. Here's a sample of their HTML documentation.

http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-A02DB70F-9704-42A4-9071-D409D783D911.htm

Solution 3

If you want to use LAPACK from C++ you might want to have a look a FLENS. It defines low- and high-level interfaces to LAPACK but also re-implements some LAPACK functions.

With the low-level FLENS-LAPACK interface you can use your own matrix/vector types (if they have a LAPACK conform memory layout). Your call of dgetrf would look like that:

info = lapack::getrf(NoTrans, dim, nrhs, a.begin(), LDA, ipiv);

and for dgetrs

lapack::getrs(NoTrans, dim, nrhs, a.begin(), LDA, ipiv, b.begin(), LDB);

So the low-level FLENS-LAPACK functions are overloaded with respect to the element types. Consequently LAPACK function sgetrs, dgetrs, cgetrs, zgetrs are in the low-level interface of FLENS-LAPACK lapack::getrs. You also pass parameters by value/reference and not as pointer (e.g. LDA instead of &LDA).

If you use the FLENS matrix-types you can code it as

info = lapack::trf(NoTrans, A, ipiv);
if (info==0) {
    lapack::trs(NoTrans, A, ipiv, b);
}

Or you just use the LAPACK driver function dgesv

lapack::sv(NoTrans, A, ipiv, b);

Here a list of FLENS-LAPACK driver functions.

Disclaimer: Yes, FLENS is my baby! That means I coded about 95% of it and every line of code was worth it.