温馨提示:本文翻译自stackoverflow.com,查看原文请点击:matlab - LAPACK dgeev function doesn't work in eigenvectors estimation

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

发布于 2020-03-30 21:35:28

我有使用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矩阵是否正确。

查看更多

提问者
jualsevi
被浏览
31
jualsevi 2020-01-31 18:04

问题是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];
    }
}

发布
问题

分享
好友

手机
浏览

扫码手机浏览