zheev 给出错误的特征值（对照 zgeev 和 numpy.linalg.eig 检查）
zheev giving wrong eigenvalues (checked against zgeev and numpy.linalg.eig)
修补 linalg 库，我试图在 c++ 中使用 LAPACKE
获得厄密矩阵向上和 运行 的对角化例程
我按照 this 示例使用 ZHEEV，然后对照其他一些方法进行检查，特别是 numpy 的 eig 和 LAPACK(E) 的 zgeev。我不想用intel自有品牌的东西，所以我避开了MKL，直接用了LAPACKE，但大部分代码和例子中的一样。
为清楚起见，我认为 general zgeev 没有理由不能处理特定情况hermitian 矩阵，即使 zheev 已优化。
这是 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" );
}
}
用
编译
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 宏，但我能找到的所有 documentation 都说我应该使用双复数，所以我有点迷茫。无论如何，如果那是问题所在，为什么 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 一致（最多有一些 rounding/machine 错误）。上面链接的英特尔教程也确认了结果。 zheev方法明明是少数，就是不知道为什么。
我已经在几台机器上试过了：
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
感谢任何帮助。
这是我在 运行 你的 python 脚本时得到的：
$ ./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 ]]
我不知道应该匹配什么。特征值匹配，但特征向量不匹配。
在编译行中将 -cblas 替换为 -blas 即可解决问题。
cblas 包一定有错误。
我刚刚不得不处理 ZGEEV() 并且偶然发现了这个问题。在这个相当旧的系统 (Ubuntu 14.04) 上，使用 BLAS 3 支持的 LAPACK，这个样本的结果也不符合预期。
所以在我的案例中，罪魁祸首是 lapack_complex_double 矩阵的初始化。我没有像原始示例那样内联，而是选择将其读入 double[2] 类型，然后使用添加的 set_matrix() 函数将其重新初始化为 lapack_complex_double。
此外，为了正确输出来自 ZGEEV 的特征向量，应将 wgr 数组传递给它，然后打印出来。
下面是更新后的源代码及其输出。希望这对探索 LAPACKE_zgeev 和 LAPACKE_zheev 函数的其他人有用（仍然没有太多文档）。
#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]);
}
}
构建和 运行：gcc sample.c -llapacke && ./a.out。输出：
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)
注意到，ZHEEV 和 ZGEEV 的特征值和特征向量确实相同，只是它们的顺序不同。因此，确保源矩阵的格式正确非常重要，因为数据最终会通过引用传递给 Fortran ('garbage-in -> garbage-out')。
修补 linalg 库，我试图在 c++ 中使用 LAPACKE
获得厄密矩阵向上和 运行 的对角化例程我按照 this 示例使用 ZHEEV，然后对照其他一些方法进行检查，特别是 numpy 的 eig 和 LAPACK(E) 的 zgeev。我不想用intel自有品牌的东西，所以我避开了MKL，直接用了LAPACKE，但大部分代码和例子中的一样。
为清楚起见，我认为 general zgeev 没有理由不能处理特定情况hermitian 矩阵，即使 zheev 已优化。
这是 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" );
}
}
用
编译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 宏，但我能找到的所有 documentation 都说我应该使用双复数，所以我有点迷茫。无论如何，如果那是问题所在，为什么 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 一致（最多有一些 rounding/machine 错误）。上面链接的英特尔教程也确认了结果。 zheev方法明明是少数，就是不知道为什么。
我已经在几台机器上试过了：
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
感谢任何帮助。
这是我在 运行 你的 python 脚本时得到的：
$ ./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 ]]
我不知道应该匹配什么。特征值匹配，但特征向量不匹配。
在编译行中将 -cblas 替换为 -blas 即可解决问题。
cblas 包一定有错误。
我刚刚不得不处理 ZGEEV() 并且偶然发现了这个问题。在这个相当旧的系统 (Ubuntu 14.04) 上，使用 BLAS 3 支持的 LAPACK，这个样本的结果也不符合预期。
所以在我的案例中，罪魁祸首是 lapack_complex_double 矩阵的初始化。我没有像原始示例那样内联，而是选择将其读入 double[2] 类型，然后使用添加的 set_matrix() 函数将其重新初始化为 lapack_complex_double。
此外，为了正确输出来自 ZGEEV 的特征向量，应将 wgr 数组传递给它，然后打印出来。
下面是更新后的源代码及其输出。希望这对探索 LAPACKE_zgeev 和 LAPACKE_zheev 函数的其他人有用（仍然没有太多文档）。
#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]);
}
}
构建和 运行：gcc sample.c -llapacke && ./a.out。输出：
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)
注意到，ZHEEV 和 ZGEEV 的特征值和特征向量确实相同，只是它们的顺序不同。因此，确保源矩阵的格式正确非常重要，因为数据最终会通过引用传递给 Fortran ('garbage-in -> garbage-out')。