problems with zhpevd from C++

[mboisso]

Hello,
I am trying to diagonalize a hermitian matrix H, such that H = conjtrans(U).W.U from C++.

I am not sure as of how to store my matrix data (column or row major), etc.

Does anyone have a simple example ? Let say I have H =
( 1, 2+3i)
(2-3i, 4)

My intuition would be to store this matrix as
{1.0,0.0,2.0,3.0,4.0,0.0}

When I call zhpevd, I get the correct eigenvalues, but I don't seem to get the correct eigenvectors.

Also, I try to do a workspace query by setting liwork, lrwork and lwork to -1. I then declare the work, iwork and rwork arrays of size 1 (as only the first element is supposed to be updated, with the optimal workspace), but I get segmentation fault. An analysis with valgrind showed me that zhpevd tries to access elements outside of these arrays.

If you have an example from C++, I would greatly appreciate.

[Michael_C_Intel4]

Hello.

For MKL LAPACK usage you should always store the matrix in column-major order. Actaully, you put down row-major order, that is, for your 2x2 testcase you should store matrix as:

{1.0,0.0,2.0,-3.0,4.0,0.0}

lwork, lrwork, liwork should be no less than minimal required by this routines - this is described in the documentation - for instance, if jobz='V' (eigenvectors are needed), lwork >= 2*n, lrwork >= 2*n^2 + 5*n + 1, liwork >= 5*n + 3. Generally they may be set larger to get more performance.
Querying workspace is a good approach to getoptimal working array sizes if you're not sure. So query a workspace, then use the corresponding first elements of working arrays to set lwork, lrwork, liwork. This is an example:

double complex zdum, *work;
double ddum, *rwork;
int idum, *iwork;
int lwork, lrwork, liwork;
...
// Workspace query, zdum used for work, ddum for rwork, idum for iwork
lwork = lrwork = liwork = -1;
zhpevd( ..., &zdum, &lwork, &ddum, &lrwork, &idum, &liwork, ... );
lwork = (int) zdum;
lrwok = (int) ddum;
liwork = idum;
work = (double complex*) malloc( lwork*sizeof(double complex) );
rwork = (double*) malloc( lrwork*sizeof(double) );
iwork = (int*) malloc( liwork*sizeof(int) );
// Make a computation
zhpevd( ..., work, &lwork, rwork, &lrwork, iwork, &liwork, ... );
...

Michael.