## 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.

## 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”

#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”