Understanding LAPACK calls in C++ with a simple example

c++, fortran, lapack

Solution

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

Problem

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.

Original source