Неточные собственные векторы из программы LAPACK ZHEEVR

Я пытался заменить восемь процедур Scipy в моем коде Cython на LAPACK уровня C. По сути, я хочу вычислить собственные векторы, соответствующие наибольшим и наименьшим собственным значениям эрмитовой матрицы размером около 50x50. Я сделал небольшой пример 3x3, чтобы убедиться, что получаю тот же результат, что и для eigh, который, кажется, внутренне использует подпрограмму ZHEEVR. Scipy предоставляет версии подпрограмм LAPACK уровня C для импорта в Cython, который я использую. Я думаю, что синтаксис и т. д. я понял правильно, поскольку один из собственных векторов и все собственные значения абсолютно одинаковы как для eigh, так и для ZHEEVR. Однако другие собственные векторы ZHEEVR неточны. На самом деле я не получаю правильного разложения.

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)

Это дает результаты

[[ 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

Инфо=0 говорит о том, что алгоритм выполнился правильно, и собственные значения верны, но Vectors_zheevr_np.T.conj()@A@vectors_zheevr_np является, как видите, не диагональной матрицей, как должно быть для корректного разложения. - верхний правый и нижний левый элементы ненулевые.

Кто-нибудь сталкивался с этой проблемой раньше, знает, как ее исправить, или есть какие-нибудь идеи? Большое спасибо!