zheev дает неправильные собственные значения (проверено на zgeev и numpy.linalg.eig)

Работая с библиотеками linalg, я попытался запустить процедуру диагонализации для эрмитовых матриц и запустить ее на c ++, используя LAPACKE.

Я последовал за этот Например, использовать ZHEEV, а затем проверить некоторые другие методы, в частности, numpy’s eig и LAPACK (E) zgeev. Я не хотел использовать материал от собственной марки Intel, поэтому я избегал MKL и просто выбрал LAPACKE, но большая часть кода такая же, как в примере.

Для ясности я не вижу причин, почему GEнерал гGEев не должен быть в состоянии справиться с конкретным случаем онматрица rmitian, даже если zонев оптимизирован.

Вот с ++

#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" );
}
}

составлено с

g++ diag.cc -L /usr/lib/lapack/ -llapacke -lcblas -o diag.out

Единственные нестандартные liblapacke-dev, а также libcblas-devустанавливается через apt-get install, Что возможно могло пойти не так?

Выход

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)

Я попытался использовать верхний треугольник, заполнить матрицу и различные другие исправления. Те же результаты каждый раз.

Я с подозрением отношусь к #define ldcmplex lapack_complex_double макрос, но все документация Я могу найти, говорит, что я должен использовать двойной комплекс, так что я немного растерялся. Во всяком случае, если бы это была проблема, почему бы работал zgeev?

Во всяком случае, вот скрипт проверки Python:

#!/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

Согласен с zgeev (до некоторых ошибок округления / машин). Результат также подтверждается указанным выше учебником Intel. Метод Жеева явно в меньшинстве одного, я просто не знаю почему.

Я пробовал это на нескольких машинах:

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

Любая помощь приветствуется.