How to start using LAPACK in c++?

Question

I'm new to computational science and I already have learned basic methods for integration, interpolation, methods like RK4, Numerov etc on c++ but recently my professor asked me to learn how to use LAPACK for solving problems related to matrices. Like for example finding eigenvalues of a complex matrix. I have never used third-party libraries and I almost always write my own functions. I have been searching around for several days but can't find any amateur-friendly guide to lapack. All of them are written in words I don't understand and I don't know why using already written functions should be this complicated. They are full of words like zgeev, dtrsv, etc. and I'm frustrated. I just want to code something like this pseudo-code:

#include <lapack:matrix>
int main(){
  LapackComplexMatrix A(n,n);
  for...
   for...
    cin>>A(i,j);
  cout<<LapackEigenValues(A);
  return 0;
}

I don't know if I'm being silly or amateur. But again, this shouldn't be that hard should it? I don't even know should I use LAPACK or LAPACK++. (I write codes in c++ and have no knowledge of Python or FORTRAN) and how to install them.

Answer 1

I'm going to disagree with some of the other answers and say that I believe that figuring out how to use LAPACK is important in the field of scientific computing.

However, there is a large learning curve to using LAPACK. This is because it is written at a very low level. The disadvantage of that is that it seems very cryptic, and not pleasant to the senses. The advantage of it is that the interface is unambiguous and basically never changes. Additionally, implementations of LAPACK, such as the Intel Math Kernel Library are really fast.

For my own purposes, I have my own higher level C++ classes which wrap around LAPACK subroutines. Many scientific libraries also use LAPACK underneath. Sometimes it's easier to just use them, but in my opinion there's a lot of value in understanding the tool underneath. To that end, I've provided a small working example written in C++ using LAPACK to get you started. This works in Ubuntu, with the liblapack3 package installed, and other necessary packages for building. It can probably be used in most Linux distributions, but installation of LAPACK and linking against it can vary.

Here's the file test_lapack.cpp

    #include <iostream>
    #include <fstream>


    using namespace std;

    // dgeev_ is a symbol in the LAPACK library files
    extern "C" {
    extern int dgeev_(char*,char*,int*,double*,int*,double*, double*, double*, int*, double*, int*, double*, int*, int*);
    }

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

      // check for an argument
      if (argc<2){
        cout << "Usage: " << argv[0] << " " << " filename" << endl;
        return -1;
      }

      int n,m;
      double *data;

      // read in a text file that contains a real matrix stored in column major format
      // but read it into row major format
      ifstream fin(argv[1]);
      if (!fin.is_open()){
        cout << "Failed to open " << argv[1] << endl;
        return -1;
      }
      fin >> n >> m;  // n is the number of rows, m the number of columns
      data = new double[n*m];
      for (int i=0;i<n;i++){
        for (int j=0;j<m;j++){
          fin >> data[j*n+i];
        }
      }
      if (fin.fail() || fin.eof()){
        cout << "Error while reading " << argv[1] << endl;
        return -1;
      }
      fin.close();

      // check that matrix is square
      if (n != m){
        cout << "Matrix is not square" <<endl;
        return -1;
      }

      // allocate data
      char Nchar='N';
      double *eigReal=new double[n];
      double *eigImag=new double[n];
      double *vl,*vr;
      int one=1;
      int lwork=6*n;
      double *work=new double[lwork];
      int info;

      // calculate eigenvalues using the DGEEV subroutine
      dgeev_(&Nchar,&Nchar,&n,data,&n,eigReal,eigImag,
            vl,&one,vr,&one,
            work,&lwork,&info);


      // check for errors
      if (info!=0){
        cout << "Error: dgeev returned error code " << info << endl;
        return -1;
      }

      // output eigenvalues to stdout
      cout << "--- Eigenvalues ---" << endl;
      for (int i=0;i<n;i++){
        cout << "( " << eigReal[i] << " , " << eigImag[i] << " )\n";
      }
      cout << endl;

      // deallocate
      delete [] data;
      delete [] eigReal;
      delete [] eigImag;
      delete [] work;


      return 0;
    }

This can be built using the command line

    g++ -o test_lapack test_lapack.cpp -llapack

This will produce an executable named test_lapack. I've set this up to read in a text input file. Here's a file named matrix.txt containing a 3x3 matrix.

3 3
-1.0 -8.0  0.0
-1.0  1.0 -5.0
 3.0  0.0  2.0

To run the program simply type

    ./test_lapack matrix.txt

at the command line, and the output should be

--- Eigenvalues ---
( 6.15484 , 0 )
( -2.07742 , 3.50095 )
( -2.07742 , -3.50095 )

---

Comments:

• You seem thrown off by the naming scheme for LAPACK. A short description is here.
• The interface for the DGEEV subroutine is here. You should be able to compare the description of the arguments there to what I've done here.
• Note the extern "C" section at the top, and that I've added an underscore to dgeev_. That's because the library was written and built in Fortran, so this is necessary to make the symbols match when linking. This is compiler and system dependent, so if you use this on Windows, it will all have to change.
• Some people might suggest using the C interface to LAPACK. They might be right, but I've always done it this way.

Answer 2

I usually resist telling people what I think they should do rather than answering their question but in this case I'm going to make an exception.

Lapack is written in FORTRAN and the API is very FORTRAN-like. There is a C API to Lapack that makes the interface slightly less painful but it will never be a pleasant experience to use Lapack from C++.

Alternatively, there is a C++ matrix class library called Eigen that has many of the capabilities of Lapack, provides computational performance comparable to the better Lapack implementations, and is very convenient to use from C++. In particular, here is how your example code might be written using Eigen

#include <iostream>
using std::cout;
using std::endl;

#include <Eigen/Eigenvalues>

int main()
{
  const int n = 4;
  Eigen::MatrixXd a(n, n); 
  a <<
    0.35, 0.45, -0.14, -0.17,
    0.09, 0.07, -0.54, 0.35,
    -0.44, -0.33, -0.03, 0.17,
    0.25, -0.32, -0.13, 0.11;
  Eigen::EigenSolver<Eigen::MatrixXd> es; 
  es.compute(a);
  Eigen::VectorXcd ev = es.eigenvalues();
  cout << ev << endl;
}

This example eigenvalue problem is a test case for the Lapack function dgeev. You can view the FORTRAN code and results for this problem dgeev example and make your own comparisons.

Answer 3

Here's another answer in the same vein as the above.

You should look into the Armadillo C++ linear algebra library.

Pros:

1. The function syntax is high-level (similar to that of MATLAB). So no DGESV mumbo-jumbo, just X = solve( A, B ) (although there is a reason behind those oddly-looking LAPACK function names...).
2. Implements various matrix decompositions (LU, QR, eigenvalues, SVD, Cholesky, etc.)
3. It is fast when used properly.
4. It is well documented.
5. Has support for sparse matrices (you will want to look into these later).
6. You can link it against your super-optimized BLAS/LAPACK libraries for optimal performance.

Here's how @BillGreene's code would look like with Armadillo:

#include <iostream>
#include <armadillo>

using namespace std;
using namespace arma;

int main()
{
   const int k = 4;
   mat A = zeros<mat>(k,k) // mat == Mat<double>

   // with the << operator...
   A <<
    0.35 << 0.45 << -0.14 << -0.17 << endr
    0.09 << 0.07 << -0.54 << 0.35  << endr
    -0.44 << -0.33 << -0.03 << 0.17 << endr
    0.25 << -0.32 << -0.13 << 0.11 << endr;

   // but using an initializer list is faster
   A = { {0.35, 0.45, -0.14, -0.17}, 
         {0.09, 0.07, -0.54, 0.35}, 
         {-0.44, -0.33, -0.03, 0.17}, 
         {0.25, -0.32, -0.13, 0.11} };

   cx_vec eigval; // eigenvalues may well be complex
   cx_mat eigvec;

   // eigenvalue decomposition for general dense matrices
   eig_gen(eigval, eigvec, A);

   std::cout << eigval << std::endl;

   return 0;
}

Answer 4

There is the SLATE library written in C++. It provides LAPACK functionality on distributed-memory systems with accelerators.