LAPACK dgeev 函数在特征向量估计中不起作用

LAPACK dgeev function doesn't work in eigenvectors estimation

我有使用 Lapack 库在 MEX 中计算特征值和特征向量的代码:

char *JOBVL="N";
char *JOBVR="V";    

mwSignedIndex LDA, LDVL, LDVR, LWORK, INFO;
LDA = dim;
LDVL = dim;
LDVR = dim;    
LWORK = 4*dim;    

double *Awork, *WR, *WI, *VL, *VR, *WORK;
Awork = (double*)mxCalloc(LDA*dim,sizeof(double));
memcpy(Awork, A, dim*dim*sizeof(double));
WR = (double *)mxCalloc(dim,sizeof(double)); 
WI = (double *)mxCalloc(dim,sizeof(double)); 
VL = (double *)mxCalloc(LDVL*dim,sizeof(double)); 
VR = (double *)mxCalloc(LDVR*dim,sizeof(double)); 
WORK = (double *)mxCalloc(LWORK,sizeof(double)); 

dgeev(JOBVL, JOBVR, &dim, Awork, &LDA, WR, WI, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);

LWORK = (mwSignedIndex)WORK[0];    
memcpy(Awork, A, dim*dim*sizeof(double));
mxFree(WORK);
WORK = (double *)mxCalloc(LWORK,sizeof(double)); 

dgeev(JOBVL, JOBVR, &dim, Awork, &LDA, WR, WI, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
LWORK = (mwSignedIndex)WORK[0];    
memcpy(Awork, A, dim*dim*sizeof(double));
mxFree(WORK);
WORK = (double *)mxCalloc(LWORK,sizeof(double)); 

dgeev(JOBVL, JOBVR, &dim, Awork, &LDA, WR, WI, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);

mxFree(Awork);
mxFree(WR);
mxFree(WI);    
mxFree(VL);
mxFree(VR);
mxFree(WORK);

memcpy(AvaW, WR, dim*sizeof(double));
memcpy(AveW, VR, dim*dim*sizeof(double));

其中AvaW和AveW分别是输出特征值和特征向量。

我的问题是:我有一个输入矩阵 A。当 A = H(H 是 NxN 维的非对称正规矩阵)时,dgeev 运行良好。我比较了 matlab 的 eig 函数和我知道正确结果的函数 f(AvaW,AveW)。 但当 A = M 时,其中 M 为:

M = [H I;
     Z Z];

I = 眼矩阵,Z 是一个零点矩阵,M 是一个 2Nx2N 维的方阵。在这种情况下,dgeev 可以很好地计算特征值,但顺序不同于 Matlab(AvaW(1:N) = AvaMatlab(N+1:2N) 和 AvaW(N+1:2N) = AvaMatlab(1:N))。使用 dgeev 计算的自动向量是错误的,与 Matlab 的自动向量不匹配并且 f(AvaW,AveW) 的结果不正确。

有人知道问题出在哪里吗?

非常感谢,问候。

PD:在所有情况下,INFO 的值为 0,这意味着解决方案是正确的。我检查 M 和 H 矩阵是否正确。

问题出在 M 处的 I 矩阵和 Z 矩阵(左下角)的位置。由于索引错误,赋值错误。不是dgeev函数的问题。

创建运行良好的 M 的代码是:

 for(i=0; i<ny; i++){
    for(j=0; j<ny; j++){
        M[2*i*ny+j] = A[i*ny+j];
        M[2*i*ny+ny+j] = Z[i*ny+j];            
        M[2*ny*ny+2*i*ny+j] = Ipos[i*ny+j];             
        M[2*ny*ny+2*i*ny+ny+j] = Z[i*ny+j];
    }
}