LAPACKE_cheev 只有 returns 特征向量的上矩阵

Question

我需要使用 LAPACKE 计算复厄尔米特矩阵的 eigenvalues/eigenvectors。我找到了函数 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 不存储整个特征向量矩阵 - 仅存储上半部分。关于上面的 link 这应该是正确的语法等

我错过了什么吗？非常感谢您的提示。

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() 旨在转置 Hermitian 矩阵，并根据 uplo 使用其上半部分或下半部分。尽管如此，cheev('V', ... ) 的输出是由输入矩阵的标准正交特征向量构成的矩阵。

其他与特征值相关的例程也会出现此问题，例如lapacke_zheev_work.c, lapacke_zheevd_work.c

例程 lapacke_zheevr_work.c and 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()。今天，该问题仅影响 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 %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);

我举报了the issue to the github repository of Lapack。

LAPACKE_cheev 只有 returns 特征向量的上矩阵

LAPACKE_cheev only returns upper matrix of eigenvectors

eigenvalue

lapack

lapacke