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.

IT++ multiplication scalar and matrix

I got error when I tried to multiply a double with a complex matrix. It turns out because IT++ 4.0.6 only implement :
operator*(const double, const cmat)

So it gives error when do : (cmat*double)

To fix this, apply the patch below:
— itpp-4.0.6/itpp/base/operators.cpp 2008-10-08 19:52:24.000000000 +0800
+++ NumericalComputation/itpp-4.0.6/itpp/base/operators.cpp 2009-03-10 14:08:55.000000000 +0800
@@ -238,6 +238,7 @@
return temp;
}

+
cmat operator*(const std::complex &s, const mat &m)
{
it_assert_debug(m.rows() > 0 && m.cols() > 0, “operator*(): Matrix of zero length”);
@@ -261,6 +262,17 @@
return temp;
}

+cmat operator*(const cmat &m, const double &s)
+{
+ it_assert_debug(m.rows() > 0 && m.cols() > 0, “operator*(): Matrix of zero length”);
+
+ cmat temp = m;
+ for (int i = 0;i < m._datasize();i++) {
+ temp._data()[i] *= (double)s;
+ }
+ return temp;
+}
+
//———————- between matrix and scalar ——————–

//———– Multiplication of a scalar and a vec —————–

and for operators.h:
— itpp-4.0.6/itpp/base/operators.h 2008-10-08 19:52:24.000000000 +0800
+++ NumericalComputation/itpp-4.0.6/itpp/base/operators.h 2009-03-10 14:40:02.000000000 +0800
@@ -654,6 +654,8 @@
*/
cmat operator/(const cmat &m, const double &s);

+cmat operator*(const cmat &m, const double &s);
+
//———————- between vec and vectors ——————–

/*!

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’

Send/receive C++ vector using MPI

I am using IT++ for some of my code and I tried to parallelize it. It came to the point where I need to sum up all the vectors from all the process. So in MPI, I will need to use MPI_Reduce. My question is since IT++ vector is not an array, can we then use MPI_Reduce with “count” parameter or do we need to define a new MPI data type to handle it?

It turns out that C++ std::vector stores the data contiguously similar to Array [1]. And IT++ vector uses std::vector of C++ standard. Hence, we can use a similar technique as Array.

this is a simple test code:

#include
#include "mpi.h"
#include

using namespace itpp;
int main(int argc, char** argv)
{
int my_rank;
int p;

MPI_Init(&argc,&argv);
MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
MPI_Comm_size(MPI_COMM_WORLD,&p);

vec test(3),result(3);
if (my_rank==0)
{
test(0)=1; test(1)=2; test(2)=3;
}
else
{
test(0)=4; test(1)=5; test(2)=6;
}
MPI_Reduce(&test(0),&result(0), test.length(), MPI_DOUBLE,
MPI_SUM, 0, MPI_COMM_WORLD);
std::cout<<"test from "<<my_rank<<" :"<<test<<"\n";
std::cout<<"result from "<<my_rank<<" :"<<result<<"\n";
MPI_Finalize();
}

The code basically creates vector “test” and store {1,2,3} in root process and {4,5,6} in other process. After that the vector is summed up in root process and displayed. To compile this, I need to include the MPI header file and IT++ header files and libs (I also used Boost library to modify my IT++ source code)

mpicxx -I/usr/lib/openmpi/include -I/home/kurniawano/local/include -I/home/kurniawano/Download/NumericalComputation/boost_1_36_0 -L/home/kurniawano/local/lib -litpp testmpivec.cpp

Then I ran the code in two processes using

mpirun -np 2 ./a.out

The output is:
test from 1 :[4 5 6]
result from 1 :[0 8.44682e+252 5.31157e+222]
test from 0 :[1 2 3]
result from 0 :[5 7 9]

So you can see that the result in root (0) is a summation of the vector from the two processes.

References:
[1] http://www.gotw.ca/publications/mill10.htm

C++ code using acos for complex number in IT++

I am currently trying out IT++ and so far it was great. It is as easy as Octave but written for C++ programming. One thing that I had trouble was that in my code I need to compute
complex acos(complex a)
However, it gives me error

for the function acos. It turns out acos for complex number was not in the C++ 98 specification, it will only be on the C++ 0x specification. So I turn to Boost library.

To modify IT++, this is what I do:
1. go to “itpp-4.0.6/itpp/base/math” folder, and modify the file “trig_hyp.h”

2. Add this line:
#include <boost/math/complex/acos.hpp>

Of course you need to have Boost library.

3. Find the line for “acos” in trig_hyp.h file, and add these lines:
//! Inverse cosine function vector
inline cvec acos(const cvec &x) { return apply_function<std::complex >(boost::math::acos, x); }
//! Inverse cosine function matrix
inline cmat acos(const cmat &x) { return apply_function<std::complex >(boost::math::acos, x); }
//! Inverse cosine function scalar
inline std::complex acos(const std::complex &x) { return boost::math::acos( x); }

4. Go back to the root directory for IT++ and do a ./configure again
./configure --prefix=$(HOME)/local F77=gfortran CPPFLAGS=-I
where we have included the root folder of our Boost libraries.

5. Type “make”, “make check”, and “make install”

6. Now, in your code, you also need to include your Boost folder in your Makefile.

That’s all, now we have acos in IT++