LAPACKE_cheev只返回特征向量的上部矩阵

Question

我需要使用 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 不存储整个特征向量矩阵 - 仅存储上半部分。关于上面的链接，这应该是正确的语法等。

我错过了什么吗？ 我将不胜感激。

Answer 1

是的，LAPACKE 的lapacke_cheev_work.c中有一个错误，在 3.9.0 版中，在行：

/* 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 );
}

英特尔软件开发工具中提供的示例也可能受到此问题的困扰，因为它也使用了LAPACKE_cheev( LAPACK_ROW_MAJOR, 'V',...)，这确实需要换位。

查看 Lapack 中的 LAPACKE 源代码，版本 3.8.0 不受此问题的影响，因为它在任何地方都使用LAPACKE_cge_trans()。今天，该问题仅影响 LAPACKE 3.9.0,2019 年 11 月 21 日。

下面是一个测试它的示例代码，它像你一样链接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 %dn",info); }   
info = LAPACKE_cheev(LAPACK_ROW_MAJOR, 'V', 'U', n, eigenvectors, lda, w);
if(info !=0){printf("LAPACKE_cheev error %dn",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 %dn",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 %dn",info); }
//need to transpose back the whole result
LAPACKE_cge_trans(LAPACK_COL_MAJOR, n, n, eigenvectors_t, n, eigenvectors,n); 

我向Lapack的github存储库报告了这个问题。