Lapack

From Eigen
Revision as of 13:43, 6 July 2010 by Ggael (Talk | contribs)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Interfacing Eigen with LAPACK

LAPACK is a large linear algebra library written in FORTRAN. It has loads of routines for all kinds of matrix problems so it is useful if you need something beyond the standard SVD, LU decomposition and so on. I needed to use it to solve the generalised eigen-problem in order to implement ellipse fitting. I shall use this as an example.

Choose your routine

All LAPACK functions have stupidly short names. Apparently it is a FORTRAN limitation. There is a list of names and what they do at the LAPACK page here: [1]. The function I wanted was dggev which "Computes the generalized eigenvalues, and left and/or right generalized eigenvectors for a pair of nonsymmetric matrices". Open the FORTRAN file for your chosen function and note the function parameters and their types. In my case it is:

     SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
    $                  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )

The types are listed in the file and generally are either INTEGER, DOUBLE PRECISION array, or CHARACTER(s).

Install LAPACK

On Ubuntu run:

sudo apt-get install liblapack

Find the Symbol

We need to check the symbol name for the function. Usually the FORTRAN compiler just adds a _ to the end. Verify like so:

$ readelf -s /usr/lib/liblapack.so | grep dggev
 1198: 0000000000228ce0 15017 FUNC    GLOBAL DEFAULT   11 dggevx_
 1350: 0000000000225e90 11852 FUNC    GLOBAL DEFAULT   11 dggev_

When you compile your program you need to link with the LAPACK library using '-llapack'.

Create the C++ Function Declaration

All FORTRAN function arguments are passed by reference (i.e. you pass a pointer to the values). So our function becomes:

extern "C" void dggev_(const char* JOBVL, const char* JOBVR, const int* N,
                       const double* A, const int* LDA, const double* B, const int* LDB,
                       double* ALPHAR, double* ALPHAI, double* BETA,
                       double* VL, const int* LDVL, double* VR, const int* LDVR,
                       double* WORK, const int* LWORK, int* INFO);

Storage order

FORTRAN matrices are stored in column-major order, just like Eigen's default matrices. In case you use row-major matrices, you can either copy them to column major matrices, or transpose them inplace:

typedef Matrix<double,Dynamic,Dynamic,RowMajor> RowMajorMatrixXd;

RowMajorMatrixXd mat;
MatrixXd col_mat = mat; // col_mat is Lapack compatible

mat.transposeInplace(); // now mat is Lapack compatible

Note also that when the LAPACK code talks about the 'order' of the matrix, they mean the number of columns.

Conclusion

That's pretty much all you need to know. The LAPACK docs tell the rest. One extra hint though: Many LAPACK functions need a scratch space to do calculations - the workspace. The size you need to make this can be found by calling the function with LWORK=-1. Here is my code for the generalised eigen-problem:


// Generalised Eigen-Problem
// Solve:
// A * v(j) = lambda(j) * B * v(j).
//
// v are the eigenvectors and are stored in v.
// lambda are the eigenvalues and are stored in lambda.
// The eigenvalues are stored as: (lambda(:, 1) + lambda(:, 2)*i)./lambda(:, 3)
//
// returns true on success.
bool GEP(const MatrixXd& A, const MatrixXd& B, MatrixXd& v, MatrixXd& lambda)
{
  int N = A.cols(); // Number of columns of A and B. Number of rows of v.
  if (B.cols() != N  || A.rows()!=N || B.rows()!=N)
    return false;

  v.resize(N,N);
  lambda.resize(N, 3);

  int LDA = A.outerStride();
  int LDB = B.outerStride();
  int LDV = v.outerStride();

  double WORKDUMMY;
  int LWORK = -1; // Request optimum work size.
  int INFO = 0;
  
  double * alphar = const_cast<double*>(lambda.col(0).data());
  double * alphai = const_cast<double*>(lambda.col(1).data());
  double * beta   = const_cast<double*>(lambda.col(2).data());

  // Get the optimum work size.
  dggev_("N", "V", &N, A.data(), &LDA, B.data(), &LDB, alphar, alphai, beta, 0, &LDV, v.data(), &LDV, &WORKDUMMY, &LWORK, &INFO);

  LWORK = int(WORKDUMMY) + 32;
  VectorXd WORK(LWORK);

  dggev_("N", "V", &N, A.data(), &LDA, B.data(), &LDB, alphar, alphai, beta, 0, &LDV, v.data(), &LDV, WORK.data(), &LWORK, &INFO);

  return INFO==0;
}