9
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Perhaps this example would be useful: matrixprogramming.com/files/code/LAPACK $\endgroup$ – nukeguy Mar 14 '17 at 13:41
  • $\begingroup$ If you are just starting, maybe it'll be easier to use a library that is simpler like ArrayFire github.com/arrayfire/arrayfire. You can directly call it from C++ and the API's are simpler and I think it can do all operations that LAPACK does. $\endgroup$ – Vikram Mar 14 '17 at 13:54
  • $\begingroup$ In this other post a user proposes his own wrapper FLENS, which has a very nice syntax that could ease your introduction to LAPACK. $\endgroup$ – Zythos Nov 1 '18 at 22:40
  • $\begingroup$ Calling LAPACK functions directly is very tedious and error prone. There are several user-friendly C++ wrappers for LAPACK which provide much easier usage, such as Armadillo. For the specific use case of complex eigen decomposition, see the user-friendly eig_gen() function, which underneath wraps this LAPACK monstrosity, zheev(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO), and reformats the obtained eigenvalues and eigenvectors into standard representations. $\endgroup$ – hbrerkere May 14 at 15:10
17
$\begingroup$

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.
$\endgroup$
  • 3
    $\begingroup$ A lot of what you're looking for can be found with some quick Googlage. Maybe you're just not sure what to search for. Netlib is the keeper of LAPACK. The documentation can be found here. This page has a handy table of the main functionality of LAPACK. Some of the important ones are (1) solving systems of equations, (2) eigenvalue problems, (3) singular value decompositions, and (4) QR factorizations. Did you understand the manual for DGEEV? $\endgroup$ – LedHead Mar 15 '17 at 12:26
  • 1
    $\begingroup$ They're all just different interfaces to the same thing. LAPACK is the original. It's written in Fortran, so to use it you have to play some games to make cross-compiling from C/C++ work, like I showed. I've never used LAPACKE, but it looks like it's a pretty thin C wrapper over LAPACK that avoids this cross compilation business, but it's still pretty low-level. LAPACK++ appears to be an even higher level C++ wrapper, but I don't think it's even supported anymore (someone correct me if I'm wrong). $\endgroup$ – LedHead Mar 15 '17 at 12:37
  • 1
    $\begingroup$ I don't know of any specific code collection. But if you Google any of the LAPACK subroutine names, you'll invariably find an old question on one of the StackExchange sites. $\endgroup$ – LedHead Mar 15 '17 at 12:41
  • 1
    $\begingroup$ @AlirezaHashemi By the way, the reason you have to provide the WORK array is because as a rule LAPACK doesn't allocate any memory inside its subroutines. If we're using LAPACK, we're likely using gobs of memory, and allocating memory is expensive, so it makes sense to let the calling routines be in charge of memory allocation. Since DGEEV requires memory to store intermediate quantities, we have to provide that working space to it. $\endgroup$ – LedHead Mar 16 '17 at 15:04
  • 1
    $\begingroup$ Got it. And I successfully wrote my first code to calculate eigenvalues of a complex matrix using zgeev. And already doing more! Thanks! $\endgroup$ – Alireza Mar 16 '17 at 15:26
7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you for your answer and explanation! I will try this library and choose the one suits my needs best. $\endgroup$ – Alireza Mar 15 '17 at 11:29
  • $\begingroup$ Oh, they overload operator,! Never seen that done in actual practice :-) $\endgroup$ – Wolfgang Bangerth Mar 16 '17 at 15:06
  • 1
    $\begingroup$ Actually, that operator, overload is more interesting/better than it might first appear. It is used to initialize matrices. The entries that initialize the matrix can be scalar constants but can also be previously-defined matrices or sub-matrices. Very MATLAB-like. Wish my C++ programming ability was good enough to implement something that sophisticated myself ;-) $\endgroup$ – Bill Greene Mar 16 '17 at 15:31
7
$\begingroup$

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;
}
$\endgroup$
  • $\begingroup$ Thank you for your answer and explanation! I will try this library and choose the one suits my needs best. $\endgroup$ – Alireza Mar 15 '17 at 11:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.