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 < itpp/base.h > 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.

Filed under: Uncategorized | Tagged: c++, C++ Libraries, IT++, numerical |

## Leave a Reply