1

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

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

4

1 回答 1

0

是的，LAPACKE 的lapacke_cheev_work.c3.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并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()无处不在。今天，该问题仅影响 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);

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

于 2019-12-30T12:52:48.577 回答