Adding Generalized Eigenvalue functions to IT++

I have added a functions call to lapack DGGEV and ZGGEV functions to solve generalized eigenvalue problem. There are three files to modify: lapack.h, eigen.h, and eigen.cpp. The diff files can be obtained from this link.

To test the code, I ran an example for generalized complex eigenvalue problem found in the NAG site. The test program is shown below:

#include 

using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

int main()
{

  cmat cA,cB;
  cA="-21.10-22.50i 53.50-50.50i -34.50+127.50i 7.50+0.50i; 
        -0.46-7.78i -3.50-37.50i -15.50+58.50i -10.50-1.50i; 
        4.30-5.50i 39.70-17.10i -68.50+12.50i -7.50-3.50i; 
        5.50+4.40i 14.40+43.30i -32.50-46.00i -19.00-32.50i";

  cB="1.00-5.00i 1.60+1.20i -3.00+0.00i 0.00-1.00i; 
        0.80-0.60i 3.00-5.00i -4.00+3.00i -2.40-3.20i;  
        1.00+0.00i  2.40+1.80i -4.00-5.00i 0.00-3.00i;  
        0.00+1.00i -1.80+2.40i  0.00-4.00i  4.00-5.00i";

  cout<< "cA = " << cA<<endl;
  cout<< "cB = " << cB<<endl;

  cvec lambda;
  cmat evecs;
  eig(cA,cB,lambda,evecs);
  cout<< "eig(cA,cB)\n";
  cout << "lambda = "<<lambda<<endl;
  cout << "evecs = "<<evecs<<endl;

  //Exit program:
  return 0;

}

And the output is:

cA = [[-21.1-22.5i 53.5-50.5i -34.5+127.5i 7.5+0.5i]
 [-0.46-7.78i -3.5-37.5i -15.5+58.5i -10.5-1.5i]
 [4.3-5.5i 39.7-17.1i -68.5+12.5i -7.5-3.5i]
 [5.5+4.4i 14.4+43.3i -32.5-46i -19-32.5i]]
cB = [[1-5i 1.6+1.2i -3+0i 0-1i]
 [0.8-0.6i 3-5i -4+3i -2.4-3.2i]
 [1+0i 2.4+1.8i -4-5i 0-3i]
 [0+1i -1.8+2.4i 0-4i 4-5i]]
eig(cA,cB)
lambda = [3-9i 2-5i 3-1i 4-5i]
evecs = [[-0.823768-0.176232i 0.639741+0.360259i 
                  0.977535+0.0224645i -0.906234+0.0937662i]
 [-0.152951+0.0706552i 0.0041597-0.000546503i 
                  0.159101-0.11371i -0.0074303+0.00687504i]
 [-0.0706552-0.152951i 0.0402123+0.0226448i 
                  0.120899-0.15371i 0.0302078-0.00312554i]
 [0.152951-0.0706552i -0.0226448+0.0402123i 
                  0.15371+0.120899i -0.0145859-0.14097i]]

Calling lapack functions from C++ codes

I have been using IT++ for my C++ class matrix and vectors. It is a great libraries. However, it lacks the function to solve generalized eigenvalue problem. so I need to link directly to lapack to do this. I follow the wrapper used in eigen.cpp to link to ZGGEV function of lapack. I tested the code with the problem shown in NAG website for ZGGEV examples.

To use IT++ matrix and vector classes, we need to include:
#include <itpp/base.h>

Then we need to declare the ZGGEV function prototype. Refer to this page for details on ZGGEV parameters. We can then define

extern "C"{
void zggev_(char *jobvl, char *jobvr, int *n, std::complex *a,
            int *lda, std::complex *b, int *ldb, std::complex *alpha,
            std::complex *beta, std::complex *vl,
            int *ldvl, std::complex *vr, int *ldvr,
            std::complex *work, int *lwork, double *rwork, int *info);
}

To call the functions, we simply type:

zggev_(&jobvl, &jobvr, &n, cA._data(), &lda, 
           cB._data(), &ldb, alpha._data(),  beta._data(), vl._data(), 
           &ldvl, vr._data(), &ldvr, work._data(), &lwork, 
           rwork._data(), &info);

and then we need to compile and link. These are the command:

gcc -I$HOME/local/include -I$HOME/Download/boost_1_36_0 \
     -L$HOME/local/lib testzlapack.cpp -o tzlapack \
     -llapack -lcblas -lf77blas -latlas -litpp -g


Note that I use a modified IT++ library that uses Boost libraries to compute complex function acos.

The complete source code is:

//file:testzlapack.cpp
#include &lt itpp/base.h &gt

using namespace itpp;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

extern "C"{
void zggev_(char *jobvl, char *jobvr, int *n, std::complex *a,
            int *lda, std::complex *b, int *ldb, std::complex *alpha,
            std::complex *beta, std::complex *vl,
            int *ldvl, std::complex *vr, int *ldvr,
            std::complex *work, int *lwork, double *rwork, int *info);
}

int main()
{

  cmat cA,cB;
  cA="-21.10-22.50i 53.50-50.50i -34.50+127.50i 7.50+0.50i; 
    -0.46-7.78i -3.50-37.50i -15.50+58.50i -10.50-1.50i; 
    4.30-5.50i 39.70-17.10i -68.50+12.50i -7.50-3.50i; 
    5.50+4.40i 14.40+43.30i -32.50-46.00i -19.00-32.50i";

 cB="1.00-5.00i 1.60+1.20i -3.00+0.00i 0.00-1.00i; 
   0.80-0.60i 3.00-5.00i -4.00+3.00i -2.40-3.20i;  
   1.00+0.00i  2.40+1.80i -4.00-5.00i 0.00-3.00i;  
   0.00+1.00i -1.80+2.40i  0.00-4.00i  4.00-5.00i";

  cout<< "cA = " << cA<<endl;
  cout<< "cB = " << cB<<endl;

  cvec lambda;
  cmat evecs;
//  eig(cA,lambda,evecs);
  char jobvl = 'N', jobvr = 'V';
  int n, lda, ldb, ldvl, ldvr, lwork, info;
  n=lda=cA.rows();
  ldb = cB.rows();
  ldvl = 1;
  ldvr = n;
  lwork = std::max(1,  n*n+64); // This may be choosen better!

  cvec work(lwork);
  vec rwork(8*n); // This may be choosen better
  cvec alpha(n), beta(n);
  cmat vl(1,1), vr(n, n);
  zggev_(&jobvl, &jobvr, &n, cA._data(), &lda, 
       cB._data(), &ldb, alpha._data(),  beta._data(), vl._data(), 
       &ldvl, vr._data(), &ldvr, work._data(), &lwork, 
       rwork._data(), &info);
  lambda=elem_div(alpha,beta);
  evecs=vr;
  cout<<endl;
  cout<< "eig(cA,cB)= \n";
  cout << "lambda = "<<lambda<<endl;
  cout << "evecs = "<<evecs<<endl;

  //Exit program:
  return 0;

}


And we get the output:

cA = [[-21.1-22.5i 53.5-50.5i -34.5+127.5i 7.5+0.5i]
 [-0.46-7.78i -3.5-37.5i -15.5+58.5i -10.5-1.5i]
 [4.3-5.5i 39.7-17.1i -68.5+12.5i -7.5-3.5i]
 [5.5+4.4i 14.4+43.3i -32.5-46i -19-32.5i]]
cB = [[1-5i 1.6+1.2i -3+0i 0-1i]
 [0.8-0.6i 3-5i -4+3i -2.4-3.2i]
 [1+0i 2.4+1.8i -4-5i 0-3i]
 [0+1i -1.8+2.4i 0-4i 4-5i]]

 eig(cA,cB)=
lambda = [3-9i 2-5i 3-1i 4-5i]
evecs = [[-0.823768-0.176232i 0.639741+0.360259i 
                     0.977535+0.0224645i -0.906234+0.0937662i]
 [-0.152951+0.0706552i 0.0041597-0.000546503i 
                     0.159101-0.11371i -0.0074303+0.00687504i]
 [-0.0706552-0.152951i 0.0402123+0.0226448i 
                     0.120899-0.15371i 0.0302078-0.00312554i]
 [0.152951-0.0706552i -0.0226448+0.0402123i 
                     0.15371+0.120899i -0.0145859-0.14097i]]

which agrees with the results shown in NAG site.

undefined reference to `vtable for

I got this error when I compile a class:
undefined reference to `vtable for Device’

it turns out that GCC requires every virtual function to be defined, so I need to add an inline {} just to satisfy this.

virtual void setGridSpc(double a){};

instead of

virtual void setGridSpc(double a);

Maintaining C++ Codes

Some Ground rules:
1. Keep class data members private.

2. About global name space:
– Avoid data with external linkage at file scope: put all global variables in a structure, class, namespace, etc. This is to avoid name collison.
– Avoid free functions at file scope in .h files, avoid free functions with external linkage in .c files.
– Avoid enumerations, typedefs, and constants at file scope in .h files.
– Avoid using preprocessor macros in .h files, except as include guards.
– Only classes, structures, unions, free operator functions should be declared at file scope in a .h file; only classes, structures, unions, and inline functions should be defined at file scope in a .h file.

3. Place a unique and predictable include guard around the contents of each header file.

4. Place a redundant include guard around each preprocessor include directive in every header file.

5. Documentation the interfaces so that they are usable by others; have at least one other developer review each interface.

6. Use a consistent method to highlight class data member, const, statics. E.g.  adding a prefix d_ to class data members, or s_ for static data, ALL CAPS for const.

7. Don’t use “using namespace” directive in header files. It forfeit the purpose of having namespaces.

taken from: Large Scale C++ Software Design by John Lakos, 1996. Thinking in C++, 2nd Edition.

Global Constants in C++

Declaring the global constants can polute the names available and create conflicts as the project grow larger. One solution is to declare it using a namespace.

globalconst.h

namespace GlobalConst{

const double PI=3.14;

};

and in your file, simply use

GlobalConst::PI

UPDATE: I have corrected the code, thanks

Compiling ARPACK++

I have had difficulties in compiling the examples code provided by ARPACK++. After some googling, it turns out that ARPACK++ is not compatible with newer version of GCC. Luckily, I found a page that provide a patch. So let me provide the steps to compile ARPACK++ .

1. First you have to install ARPACK, UMFPACK, and SuperLU.

2. To install ARPACK, download the files from (http://www.caam.rice.edu/software/ARPACK/) make sure to download the “patch” file as well. Then simply extract the two files:

        zcat arpack96.tar.gz | tar -xvf -
        zcat patch.tar.gz    | tar -xvf -
3. In the ARPACK folder, modify ARmake.inc. It is better to use the BLAS and LAPACK provided by ARPACK. Sometimes using other libraries creates problem. And for my case, I need to modify the location of the ARPACK source tree and the location for "make" program.
4. type make lib. This should compile arpack for your platform, copy the library to your local library locations. You should test also the EXAMPLES codes provided to make sure everything went well.
5. Now install SuperLU, downlaod the source code from (http://crd.lbl.gov/~xiaoye/SuperLU/#superlu).
6. Modify make.inc. This is the one that I use:
SuperLUroot     = /scratch/kurniawano/Download/SuperLU_3.1
SUPERLULIB      = $(SuperLUroot)/lib/libsuperlu_3.1.a
BLASDEF         = -DUSE_VENDOR_BLAS
BLASLIB         = /usr/lib64/libblas.a
TMGLIB          = libtmglib.a
LIBS            = $(SUPERLULIB) $(BLASLIB)
and
CDEFS        = -DAdd_
with "one" underscore after "-DAdd"
7. Then type make, if there is an error due to the timer, just type make again and it will continue compiling.
8. Now UMFPACK is installed from the SuiteSparse package. I have written some other posts on this package installation. Just search on this blog.
9. So now download Arpack++ from (http://www.ime.unicamp.br/~chico/arpack++/), and its patch from (http://reuter.mit.edu/index.php/software/arpackpatch/).
10. extract Arpack++ and copy the patch file to the arpack++ folder, and type 
patch -p1 < arpack++1.2.patch.diff
11. After that modify Makefile.inc, especially the path to the Arpack++ and the following libraries:

ARPACK_LIB   = /scratch/kurniawano/local/lib/libarpack.a
LAPACK_LIB   = -llapack
UMFPACK_LIB  = /scratch/kurniawano/local/lib/libumfpack.a
SUPERLU_LIB  = /scratch/kurniawano/local/lib/libsuperlu.a
BLAS_LIB     = -lblas
FORTRAN_LIBS = /scratch/kurniawano/local/lib/libg2c.a

Compiling IT++ with MKL library

I downloaded it++ 4.0.6 and tried to compile using MKL library. To do that, I need first to set the environment path

export LDFLAGS=-L/opt/intel/Compiler/11.0/081/mkl/lib/32/

export CPPFLAGS=-I/opt/intel/Compiler/11.0/081/mkl/include

and then I run the configure with this option

./configure –prefix=$HOME/local F77=gfortran –with-blas=”-lmkl_intel -lmkl_sequential -lmkl_core -lpthread” –with-lapack=”-lmkl_lapack”

With this settings configure can detect my MKL library

UPDATE: though it can detect MKL, but it gives a compilation error. I haven’t managed to solve it. so I currently by passing the mkl library and install ATLAS instead.

./configure –prefix=/scratch/ihpcoka/local CPPFLAGS=’-I$(HOME)/DATA/Download/boost_1_36_0 -I/scratch/ihpcoka/local/include’ F77=gfortran LDFLAGS=’-L/usr/lib64 -L/scratch/ihpcoka/local/lib’ –with-fftw=’-lfftw3′ –with-blas=’-lf77blas -latlas’ –with-lapack=’-llapack’