Reputation: 31
zheev giving wrong eigenvalues (checked against zgeev and numpy.linalg.eig)
Tinkering with linalg libraries, I tried to get a diagonalisation routine for hermitian matrices up and running in c++ using LAPACKE
I followed this example to use ZHEEV, and then checked against some other methods, specifically numpy's eig and LAPACK(E)'s zgeev. I didn't want to use intel's own-brand stuff so I avoided MKL and just went for straight LAPACKE, but the majority of the code is the same as in the example.
For clarity, I see no reason why the general zgeev shouldn't be able to handle the specific case of a hermitian matrix, even if zheev is optimised.
Here's the c++
#include <stdlib.h>
#include <stdio.h>
#include <lapacke.h>
//Parameters
#define N 4
#define LDA N
#define lint lapack_int
#define ldcmplex lapack_complex_double
//Auxiliary routines prototypes
extern void print_matrix( char* desc, lint m, lint n, ldcmplex* a, lint lda );
extern void print_rmatrix( char* desc, lint m, lint n, double* a, lint lda );
//Main program
int main()
{
//Locals
lint n = N, lda = LDA, info;
;
//Local arrays
double wr[N];
ldcmplex ah[LDA*N] = {
{ 9.14, 0.00}, { 0.00, 0.00}, { 0.00, 0.00}, { 0.00, 0.00},
{-4.37, 9.22}, {-3.35, 0.00}, { 0.00, 0.00}, { 0.00, 0.00},
{-1.98, 1.72}, { 2.25, 9.51}, {-4.82, 0.00}, { 0.00, 0.00},
{-8.96, 9.50}, { 2.57, -2.40}, {-3.24, -2.04}, { 8.44, 0.00}
};
;
//Executable statements
printf( "LAPACKE_zheev (row-major, high-level) Example Program Results\n" ) ;
;
//Print martix
print_matrix( "Input Matrix", n, n, ah, lda );
;
//Solve eigenproblem
info = LAPACKE_zheev( LAPACK_ROW_MAJOR, 'V', 'L', n, ah, lda, wr );
;
//Check for convergence
if( info > 0 ) {
printf( "zheev algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
;
//Print eigenvalues
print_rmatrix( "zheev Eigenvalues", 1, n, wr, 1 );
;
//Print eigenvectors
print_matrix( "Eigenvectors (stored columnwise)", n, n, ah, lda );
;
//Local arrays
ldcmplex wc[N];
ldcmplex ag[LDA*N] = {
{ 9.14, 0.00}, {-4.37, -9.22}, {-1.98, -1.72}, {-8.96, -9.50},
{-4.37, 9.22}, {-3.35, 0.00}, { 2.25, -9.51}, { 2.57, 2.40},
{-1.98, 1.72}, { 2.25, 9.51}, {-4.82, 0.00}, {-3.24, 2.04},
{-8.96, 9.50}, { 2.57, -2.40}, {-3.24, -2.04}, { 8.44, 0.00},
};
;
//Executable statements
printf( "LAPACKE_zgeev (row-major, high-level) Example Program Results\n" );
;
//Print martix
print_matrix( "Input Matrix", n, n, ag, lda );
;
//Solve eigenproblem
info = LAPACKE_zgeev( LAPACK_ROW_MAJOR, 'N', 'V', n, ag, lda, wc, 0, lda, 0, lda);
;
//Check for convergence
if( info > 0 ) {
printf( "zgeev algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
;
//Print eigenvalues
print_matrix( "zgeev Eigenvalues", 1, n, wc, 1);
;
//Print eigenvectors
print_matrix( "Eigenvectors (stored columnwise)", n, n, ag, lda );
exit( 0 );
}
//Auxiliary routine: printing a matrix
void print_matrix( char* desc, lint m, lint n, ldcmplex* a, lint lda ) {
lint i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ )
printf( " (%6.2f,%6.2f)", creal(a[i*lda+j]), cimag(a[i*lda+j]) );
printf( "\n" );
}
}
//Auxiliary routine: printing a real matrix
void print_rmatrix( char* desc, lint m, lint n, double* a, lint lda ) {
lint i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %6.2f", a[i*lda+j] );
printf( "\n" );
}
}
compiled with
g++ diag.cc -L /usr/lib/lapack/ -llapacke -lcblas -o diag.out
The only nonstandard packages being liblapacke-dev, and libcblas-dev, installed via apt-get install. What could possibly go wrong?
The output is
LAPACKE_zheev (row-major, high-level) Example Program Results
Input Matrix
( 9.14, 0.00) ( 0.00, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( -4.37, 9.22) ( -3.35, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( -1.98, 1.72) ( 2.25, 9.51) ( -4.82, 0.00) ( 0.00, 0.00)
( -8.96, 9.50) ( 2.57, -2.40) ( -3.24, -2.04) ( 8.44, 0.00)
zheev Eigenvalues
-18.96 -12.85 18.78 30.71
Eigenvectors (stored columnwise)
( 0.16, 0.00) ( 0.57, 0.00) ( -0.73, 0.00) ( 0.35, 0.00)
( 0.26, -0.81) ( 0.17, -0.25) ( 0.22, -0.38) ( 0.06, -0.02)
( 0.29, 0.27) ( -0.11, -0.30) ( -0.26, -0.42) ( -0.50, -0.50)
( -0.21, 0.23) ( 0.50, -0.49) ( 0.18, -0.09) ( -0.33, 0.51)
LAPACKE_zgeev (row-major, high-level) Example Program Results
Input Matrix
( 9.14, 0.00) ( -4.37, -9.22) ( -1.98, -1.72) ( -8.96, -9.50)
( -4.37, 9.22) ( -3.35, 0.00) ( 2.25, -9.51) ( 2.57, 2.40)
( -1.98, 1.72) ( 2.25, 9.51) ( -4.82, 0.00) ( -3.24, 2.04)
( -8.96, 9.50) ( 2.57, -2.40) ( -3.24, -2.04) ( 8.44, 0.00)
zgeev Eigenvalues
( 25.51, 0.00) (-16.00, -0.00) ( -6.76, 0.00) ( 6.67, 0.00)
Eigenvectors (stored columnwise)
( 25.51, 0.00) ( -0.00, 0.00) ( 0.00, 0.00) ( 0.00, -0.00)
( 0.00, 0.00) (-16.00, -0.00) ( 0.00, -0.00) ( 0.00, -0.00)
( 0.00, 0.00) ( 0.00, 0.00) ( -6.76, 0.00) ( -0.00, -0.00)
( 0.00, 0.00) ( 0.00, 0.00) ( 0.00, 0.00) ( 6.67, 0.00)
I tried using the Upper triangle, filling the matrix and various other fixes. Same results each time.
I'm suspicious of the #define ldcmplex lapack_complex_double macro, but all the documentation I can find says I should be using a double complex, so I'm a bit lost. Anyway, if that was the problem, why would zgeev work?
Anyway, here's the python check script:
#!/usr/bin/env python
from numpy import linalg as li
import numpy as np
mat=np.array([
[ 9.14 + 0.00j, 0.00 + 0.00j, 0.00 + 0.00j, 0.00 +0.00j],
[ -4.37 + 9.22j, -3.35 + 0.00j, 0.00 + 0.00j, 0.00 +0.00j],
[ -1.98 + 1.72j, 2.25 + 9.51j, -4.82 + 0.00j, 0.00 +0.00j],
[ -8.96 + 9.50j, 2.57 - 2.40j, -3.24 - 2.04j, 8.44 +0.00j]])
mat[0]=np.conj(mat[:,0])
mat[1]=np.conj(mat[:,1])
mat[2]=np.conj(mat[:,2])
mat[3]=np.conj(mat[:,3])
mat=np.matrix(mat)
w, v = li.eig(mat)
print w
print v
It agrees with zgeev (up to some rounding/machine errors). The result is also confirmed by the intel tutorial linked above. The zheev method is clearly in a minority of one, I just don't know why.
I've tried this on a couple of machines:
Linux parabox 4.8.0-52-generic #55-Ubuntu SMP Fri Apr 28 13:28:50 UTC 2017 x86_64 x86_64 x86_64 GNU/Linux
Linux glass 4.10.0-21-generic #23-Ubuntu SMP Fri Apr 28 16:14:22 UTC 2017 x86_64 x86_64 x86_64 GNU/Linux
Any and all help appreciated.
Upvotes: 1
Views: 1377
Answers (3)
Reputation: 1
I just had to deal with ZGEEV() and stumbled on this issue. On this rather old system (Ubuntu 14.04), with BLAS 3 backed LAPACK, the results from this sample were also not matching to the expected.
So it turned out the culprit in my case was the initialization of the lapack_complex_double matrix. Instead of doing it inline, as is in the original sample, I opted for reading it into a double[2] type and then reinitializing it into lapack_complex_double using the added set_matrix() function.
Additionally, for proper output of the eigen vectors from ZGEEV, a wgr array should be passed to it, and subsequently printed out.
Below is the updated source code and its output. Hope this would be useful for someone else exploring the LAPACKE_zgeev and LAPACKE_zheev functions (not much documentation still).
#include <stdlib.h>
#include <stdio.h>
#include <lapacke.h>
//Parameters
#define N 4
#define LDA N
#define lint lapack_int
#define ldcmplex lapack_complex_double
typedef struct double2 {
double v[2];
} double2_t;
//Auxiliary routines prototypes
extern void print_matrix( char* desc, lint m, lint n, ldcmplex* a, lint lda);
extern void print_rmatrix( char* desc, lint m, lint n, double* a, lint lda);
extern void set_matrix(lapack_int n, lapack_complex_double* a, lapack_int lda, double2_t *a2);
//Main program
int main()
{
//Locals
lint n = N, lda = LDA, info;
;
//Local arrays
double wr[N];
ldcmplex ah[LDA*N] = {0};
double2_t ah2[LDA*N] = {
{ 9.14, 0.00}, { 0.00, 0.00}, { 0.00, 0.00}, { 0.00, 0.00},
{-4.37, 9.22}, {-3.35, 0.00}, { 0.00, 0.00}, { 0.00, 0.00},
{-1.98, 1.72}, { 2.25, 9.51}, {-4.82, 0.00}, { 0.00, 0.00},
{-8.96, 9.50}, { 2.57, -2.40}, {-3.24, -2.04}, { 8.44, 0.00}
};
;
//Executable statements
set_matrix(n, ah, lda, ah2);
printf( "LAPACKE_zheev (row-major, high-level) Example Program Results\n" ) ;
;
//Print martix
print_matrix( "Input Matrix", n, n, ah, lda );
;
//Solve eigenproblem
info = LAPACKE_zheev( LAPACK_ROW_MAJOR, 'V', 'L', n, ah, lda, wr );
;
//Check for convergence
if( info > 0 ) {
printf( "zheev algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
;
//Print eigenvalues
print_rmatrix( "zheev Eigenvalues", 1, n, wr, 1 );
;
//Print eigenvectors
print_matrix( "Eigenvectors (stored columnwise)", n, n, ah, lda );
;
//Local arrays
ldcmplex wc[N];
ldcmplex ag[LDA*N] = {0};
ldcmplex wgr[LDA*N] = {0};
double2_t ag2[LDA*N] = {
{ 9.14, 0.00}, {-4.37, -9.22}, {-1.98, -1.72}, {-8.96, -9.50},
{-4.37, 9.22}, {-3.35, 0.00}, { 2.25, -9.51}, { 2.57, 2.40},
{-1.98, 1.72}, { 2.25, 9.51}, {-4.82, 0.00}, {-3.24, 2.04},
{-8.96, 9.50}, { 2.57, -2.40}, {-3.24, -2.04}, { 8.44, 0.00},
};
;
printf("\n\n");
//Executable statements
set_matrix(n, ag, lda, ag2);
printf( "LAPACKE_zgeev (row-major, high-level) Example Program Results\n" );
;
//Print martix
print_matrix( "Input Matrix", n, n, ag, lda );
;
//Solve eigenproblem
info = LAPACKE_zgeev( LAPACK_ROW_MAJOR, 'N', 'V', n, ag, lda, wc, 0, lda, wgr, lda);
;
//Check for convergence
if( info > 0 ) {
printf( "zgeev algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
;
//Print eigenvalues
print_matrix( "zgeev Eigenvalues", 1, n, wc, 1);
;
//Print eigenvectors
print_matrix( "Eigenvectors (stored columnwise)", n, n, wgr, lda );
exit( 0 );
}
//Auxiliary routine: printing a matrix
void print_matrix( char* desc, lint m, lint n, ldcmplex* a, lint lda ) {
lint i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ )
printf( " (%6.2f,%6.2f)", creal(a[i*lda+j]), cimag(a[i*lda+j]) );
printf( "\n" );
}
}
//Auxiliary routine: printing a real matrix
void print_rmatrix( char* desc, lint m, lint n, double* a, lint lda ) {
lint i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %6.2f", a[i*lda+j] );
printf( "\n" );
}
}
//Auxiliary routine: set a complex matrix from a double[2] type matrix
void set_matrix(lapack_int n, lapack_complex_double* a, lapack_int lda, double2_t *a2) {
lapack_int i, j;
for( i = 0; i < n; i++ ) {
for( j = 0; j < n; j++ )
a[i*lda+j] = lapack_make_complex_double(a2[i*lda+j].v[0], a2[i*lda+j].v[1]);
}
}
Build and run: gcc sample.c -llapacke && ./a.out. The output:
LAPACKE_zheev (row-major, high-level) Example Program Results
Input Matrix
( 9.14, 0.00) ( 0.00, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( -4.37, 9.22) ( -3.35, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( -1.98, 1.72) ( 2.25, 9.51) ( -4.82, 0.00) ( 0.00, 0.00)
( -8.96, 9.50) ( 2.57, -2.40) ( -3.24, -2.04) ( 8.44, 0.00)
zheev Eigenvalues
-16.00 -6.76 6.67 25.51
Eigenvectors (stored columnwise)
( 0.34, -0.00) ( -0.55, 0.00) ( 0.31, 0.00) ( -0.70, 0.00)
( 0.44, -0.54) ( 0.26, 0.18) ( 0.45, 0.29) ( 0.22, -0.28)
( -0.48, -0.37) ( -0.52, -0.02) ( -0.05, 0.57) ( 0.15, 0.08)
( 0.10, -0.12) ( -0.50, 0.28) ( -0.23, -0.48) ( 0.34, -0.49)
LAPACKE_zgeev (row-major, high-level) Example Program Results
Input Matrix
( 9.14, 0.00) ( -4.37, -9.22) ( -1.98, -1.72) ( -8.96, -9.50)
( -4.37, 9.22) ( -3.35, 0.00) ( 2.25, -9.51) ( 2.57, 2.40)
( -1.98, 1.72) ( 2.25, 9.51) ( -4.82, 0.00) ( -3.24, 2.04)
( -8.96, 9.50) ( 2.57, -2.40) ( -3.24, -2.04) ( 8.44, 0.00)
zgeev Eigenvalues
( 25.51, -0.00) (-16.00, -0.00) ( -6.76, -0.00) ( 6.67, -0.00)
Eigenvectors (stored columnwise)
( 0.70, 0.00) ( 0.22, 0.27) ( 0.48, 0.26) ( -0.03, -0.31)
( -0.22, 0.28) ( 0.70, 0.00) ( -0.14, -0.29) ( 0.24, -0.48)
( -0.15, -0.08) ( -0.01, -0.61) ( 0.45, 0.27) ( 0.58, 0.00)
( -0.34, 0.49) ( 0.16, -0.00) ( 0.58, 0.00) ( -0.46, 0.27)
As one notices, the eigenvalues and eigen vectors from both ZHEEV and ZGEEV are indeed the same, only their order is different. So it's importatnt to make sure the source matrix is properly formed, as the data eventually gets passed to Fortran by reference ('garbage-in -> garbage-out').
Upvotes: 0
Reputation: 31
Replacing -cblas with -blas in the compile line solves the issue.
The cblas package must have a bug.
Upvotes: 1
Reputation: 6158
This is what I get when I run your python script:
$ ./diag.py
[ 25.51400517 +1.20330583e-15j -16.00474647 -2.91871119e-15j
-6.76497015 -6.59630730e-16j 6.66571145 +1.54590036e-16j]
[[ 0.69747489+0.j 0.21857873+0.26662122j 0.47736933+0.26449375j -0.02829581-0.30988089j]
[-0.21578745+0.28003172j 0.69688890+0.j -0.14143627-0.2852389j 0.24437193-0.47778739j]
[-0.14607303-0.08302697j -0.01445974-0.60818924j 0.44678181+0.26546077j 0.57583205+0.j ]
[-0.34133591+0.49376693j 0.15930699-0.00061647j 0.57507627+0.j -0.45823952+0.2713093j ]]
I don't know what that's supposed to match. The eigenvalues match, but not the eigenvectors.
Upvotes: 0
Related Questions
- Lapack, DGEEV does not give me the right eigenvalues
- Incorrect eigenvalues with ZGEEV in LAPACK when repeated
- lapack zheevd gives wrong results
- Inaccurate zheev eigen values and vectors
- Incorrect/inconsistent results from zgeev() LAPACK
- Different eigenvectors are obtained for lower and upper input in zheev in lapack
- Incorrect eigenvalues with cgeev
- LAPACK giving me incorrect eigenvalues
- LAPACK routine ZGEEV - gives wrong eigenvalues
- FORTRAN ZGEEV, all 0 eigenvalues