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.
Replacing -cblas with -blas in the compile line solves the issue.
The cblas package must have a bug.