LAPACKE_cheev仅返回特征向量的上矩阵
我需要使用LAPACKE计算复数厄米特矩阵的特征值/特征向量。我找到了函数LAPACKE_cheev。它可以正确地计算特征值。然而，它只存储特征向量的上矩阵。我遵循了[https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/lapacke_cheev_row.c.htm]上的示例代码
我的代码基本上是一样的：
lapack_complex_float *eigenvectors = (lapack_complex_float*) malloc(num_receivers*num_receivers* sizeof(lapack_complex_float));
//copies upper matrix 'R' into complex matrix 'eigenvalues'
info= LAPACKE_clacpy(LAPACK_ROW_MAJOR,'U', num_receivers,num_receivers,R,num_receivers,eigenvectors,num_receivers);
int n = num_receivers;
int lda =n;
float w[n];
info = LAPACKE_cheev(LAPACK_ROW_MAJOR, 'V', 'U', n, eigenvectors, lda, w);矩阵特征向量只存储矩阵R的上半部分--这工作得很好。然而，cheev并不存储整个特征向量矩阵--只存储上半部分。关于上面的链接，这应该是正确的语法等。
我是不是遗漏了什么？如果您能给我一个提示，我将非常感激。
回答 1
Stack Overflow用户
发布于 2019-12-30 20:52:49
是的，在LAPACKE的lapacke_cheev_work.c中，3.9.0版的第行有一个bug：
/* Transpose output matrices */
LAPACKE_che_trans( LAPACK_COL_MAJOR, uplo, n, a_t, lda_t, a, lda );实际上，LAPACKE_che_trans()被设计为转置厄米特矩阵，并根据uplo利用它们的上半部分或下半部分。然而，cheev('V', ... )的输出是由输入矩阵的正交特征向量组成的矩阵。
这个问题也出现在其他与特征值相关的例程中，例如lapacke_zheev_work.c、lapacke_zheevd_work.c
例程lapacke_zheevr_work.c和lapacke_zheevx_work.c正确地使用了LAPACKE_zge_trans()
if( LAPACKE_lsame( jobz, 'v' ) ) {
LAPACKE_zge_trans( LAPACK_COL_MAJOR, n, ncols_z, z_t, ldz_t, z,ldz );
}Intel software development tools提供的示例也可能受到此问题的困扰，因为它还使用了LAPACKE_cheev( LAPACK_ROW_MAJOR, 'V',...)，这确实需要转置。
看看Lapack中的LAPACKE的源代码，3.8.0版没有受到这个问题的影响，因为它到处都使用LAPACKE_cge_trans()。今天，这个问题只影响2019年11月21日的LAPACKE 3.9.0。
下面是测试它的示例代码，您可以将LAPACKE_clacpy(LAPACK_ROW_MAJOR,'U',....)和LAPACKE_cheev(LAPACK_ROW_MAJOR, 'V', 'U',...)链接起来：
#include <stdio.h>
#include <complex.h>
#include <lapacke.h>
int main()
{
lapack_int i,j, n, lda, info;
n=4;
lda=n;
float w[n];
float complex R [n*n];
for (i=0;i<n;i++){
for (j=0;j<n;j++){
R[i*n+j]=0.;
}
R[i*n+i]=2.;
if(i>0){
R[i*n+(i-1)]=-1;
}
if(i<n-1){
R[i*n+(i+1)]=-1;
}
}
float complex eigenvectors [n*n];
info= LAPACKE_clacpy(LAPACK_ROW_MAJOR,'U', n,n,R,n,eigenvectors,n);
if(info !=0){printf("LAPACKE_clacpy error %d\n",info); }
info = LAPACKE_cheev(LAPACK_ROW_MAJOR, 'V', 'U', n, eigenvectors, lda, w);
if(info !=0){printf("LAPACKE_cheev error %d\n",info); }
for (i=0;i<n;i++){
for (j=0;j<n;j++){
printf("%+6.4f+I*%+6.4f | ",creal(eigenvectors[i*n+j]),cimag(eigenvectors[i*n+j]));
}
printf("\n");
}
if(creal(eigenvectors[(n-1)*n+0])==0.){
printf("LAPACKE_cheev : the eigenvectors are wrong due to LAPACKE_che_trans()\n");
return 1;
}
return 0;
}它是通过使用
gcc main11.c -o main11b -L/home/xxxx/softs/lapack-3.9.0/lapack-3.9.0 -llapacke -llapack -lblas -lm -lgfortran -Wall如果删除-L/home/xxxx/softs/lapack-3.9.0/lapack-3.9.0，从而使用低于3.9.0的LAPACK版本，它将输出：
-0.3717+I*+0.0000 | +0.6015+I*+0.0000 | +0.6015+I*+0.0000 | -0.3717+I*+0.0000 |
-0.6015+I*+0.0000 | +0.3717+I*+0.0000 | -0.3717+I*+0.0000 | +0.6015+I*+0.0000 |
-0.6015+I*+0.0000 | -0.3717+I*+0.0000 | -0.3717+I*+0.0000 | -0.6015+I*+0.0000 |
-0.3717+I*+0.0000 | -0.6015+I*+0.0000 | +0.6015+I*+0.0000 | +0.3717+I*+0.0000 |如果出现问题，它将打印：
-0.3717+I*+0.0000 | +0.6015+I*+0.0000 | +0.6015+I*+0.0000 | -0.3717+I*+0.0000 |
+0.0000+I*+0.0000 | +0.3717+I*+0.0000 | -0.3717+I*+0.0000 | +0.6015+I*+0.0000 |
+0.0000+I*+0.0000 | +0.0000+I*+0.0000 | -0.3717+I*+0.0000 | -0.6015+I*+0.0000 |
+0.0000+I*+0.0000 | +0.0000+I*+0.0000 | +0.0000+I*+0.0000 | +0.3717+I*+0.0000 |
LAPACKE_cheev : the eigenvectors are wrong due to LAPACKE_che_trans()可以通过重现在以前版本的lapacke_cheev_work.c中执行的步骤来纠正该问题
float complex eigenvectors_t [n*n];
info= LAPACKE_clacpy(LAPACK_ROW_MAJOR,'U', n,n,R,n,eigenvectors,n);
if(info !=0){printf("LAPACKE_clacpy error %d\n",info); }
LAPACKE_cge_trans(LAPACK_ROW_MAJOR, n, n, eigenvectors, n, eigenvectors_t,n);
info = LAPACKE_cheev(LAPACK_COL_MAJOR, 'V', 'U', n, eigenvectors_t, lda, w);
if(info !=0){printf("LAPACKE_cheev error %d\n",info); }
//need to transpose back the whole result
LAPACKE_cge_trans(LAPACK_COL_MAJOR, n, n, eigenvectors_t, n, eigenvectors,n); https://stackoverflow.com/questions/59520534
复制相似问题
腾讯云开发者
