I was trying to replace Scipy's eigh routine in my cython code with C-level LAPACK. Essentially I want to compute eigenvectors corresponding to the largest and smallest eigenvalues of a hermitian matrix of dimension around 50x50. I made a small 3x3 example to make sure I get the same result as from eigh, which internally uses the ZHEEVR routine it seems. Scipy provides C-level versions of the LAPACK routines to import in cython, which I am using. I think I got the syntax and so on right, since one of the eigenvectors and all the eigenvalues are exactly the same for both eigh and ZHEEVR. However, the other eigenvectors from ZHEEVR are inaccurate. In fact, I don't obtain a correct decomposition.
import numpy as np
cimport numpy as cnp
from scipy.linalg import eigh
from scipy.linalg.cython_lapack cimport zheevr
from libc.stdlib cimport malloc, free
dim = 3
A = (np.arange(dim**2) + 2j).reshape(dim, dim) + 2*np.eye(dim)
A = A @ np.transpose(np.conjugate(A))
print(A)
print('')
values_eigh, vectors_eigh = eigh(A)
print(values_eigh)
print('')
print(vectors_eigh)
print('')
cdef:
char jobz = 'V'
char rrange = 'A'
char uplo = 'U'
int n = dim
double complex * A_pointer = <double complex *>cnp.PyArray_DATA(
A.astype(np.complex128))
int lda = dim
double vl = 0
double vu = 0
int il = 0
int iu = 0
double abstol = 1e-12
int m = n
double * values_zheever = <double *> malloc(
n * sizeof(double))
double complex * vectors_zheever = <double complex *>malloc(
n * n * sizeof(double complex))
int ldz = n
int * isuppz = <int *> malloc(2 * n * sizeof(int))
double complex * work = <double complex *>malloc(
1 * sizeof(double complex))
int lwork = -1
double * rwork = <double *>malloc(1 * sizeof(double))
int lrwork = -1
int * iwork = <int *>malloc(1 * sizeof(int))
int liwork = -1
int info
zheevr(&jobz, &rrange, &uplo, &n, A_pointer, &lda, &vl, &vu, &il, &iu, &abstol, &m, values_zheever, vectors_zheever, &ldz, isuppz, work, &lwork, rwork, &lrwork, iwork, &liwork, &info)
lwork = <int> work[0]
lrwork = <int> rwork[0]
liwork = <int> iwork[0]
free(work)
free(rwork)
free(iwork)
work = <double complex *>malloc(lwork * sizeof(double complex))
rwork = <double *>malloc(lrwork * sizeof(double))
iwork = <int *>malloc(liwork * sizeof(int))
zheevr(&jobz, &rrange, &uplo, &n, A_pointer, &lda, &vl, &vu, &il, &iu, &abstol, &m, values_zheever, vectors_zheever, &ldz, isuppz, work, &lwork, rwork, &lrwork, iwork, &liwork, &info)
for i in range(dim):
print(values_zheever[i])
vectors_zheevr_np = np.zeros((n* n), dtype=np.complex128)
for i in range(n*n):
vectors_zheevr_np[i] = vectors_zheever[i]
vectors_zheevr_np = vectors_zheevr_np.reshape(n, n).T
print('')
print(vectors_zheevr_np)
print('')
print(vectors_zheevr_np.T.conj() @ A @ vectors_zheevr_np)
print('')
print(info)
This gives the results
[[ 21. +0.j 34.+18.j 51.+36.j]
[ 34.-18.j 82. +0.j 122.+18.j]
[ 51.-36.j 122.-18.j 197. +0.j]]
[ 0.82663284 4. 295.17336716]
[[ 0.8755536 -0.00000000e+00j -0.40824829-0.00000000e+00j
-0.25833936-0.00000000e+00j]
[ 0.24467905+6.98442317e-02j 0.81649658+1.22124533e-15j
-0.46103588+2.36713326e-01j]
[-0.38619551+1.39688463e-01j -0.40824829-6.10622664e-16j
-0.6637324 +4.73426651e-01j]]
0.8266328441181744
4.000000000000024
295.1733671558813
[[-0.8233493 +2.97808740e-01j -0.40824829-2.22044605e-16j
0.21031944-1.50016524e-01j]
[-0.20633357+1.48904370e-01j 0.81649658-1.56819002e-15j
0.51279727-7.50082620e-02j]
[ 0.41068216-0.00000000e+00j -0.40824829-0.00000000e+00j
0.8152751 +0.00000000e+00j]]
[[ 2.72444561e+01+3.55271368e-15j -1.35003120e-13-1.42108547e-14j
1.95681896e+01-8.18241075e+01j]
[-1.30784272e-13+2.16125482e-14j 4.00000000e+00-4.07921987e-16j
-2.66453526e-14+4.85653027e-13j]
[ 1.95681896e+01+8.18241075e+01j -2.84217094e-14-4.84945417e-13j
2.68755544e+02+0.00000000e+00j]]
0
The info = 0 tells me the algorithm was executed correctly, and the eigenvalues are correct, but vectors_zheevr_np.T.conj() @ A @ vectors_zheevr_np is, as you can see, not a diagonal matrix, as it should be for a correct decomposition - the upper right and bottom left elements are non-zero.
Did anyone encounter this issue before, knows how to fix it, or has any ideas? Thanks a lot!
As Nick ODell pointed out, LAPACK requires FORTRAN-contiguous input arrays. So, replacing
cdef double complex * A_pointer = <double complex *>cnp.PyArray_DATA(
A.astype(np.complex128))
with
cdef double complex * A_pointer = <double complex *>cnp.PyArray_DATA(
np.asfortranarray(A.astype(np.complex128)))
solves the problem. Furthermore, note that eigenvectors of hermitian matrices, in the cases of distinct eigenvalues, are only unique up to a multiplicative (complex) constant.