LAPACKE_cheev only returns upper matrix of eigenvectors

I need to compute the eigenvalues/eigenvectors of a complex hermitian matrix by using LAPACKE. I found the function LAPACKE_cheev. It computes the eigenvalues correctly. However, it only stores the upper matrix of eigenvectors. I followed the example code found on: [https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/lapacke_cheev_row.c.htm]

My code is basically the same:

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);

The matrix eigenvectors simply stores the upper half of the matrix R - this works just fine. However, cheev does not store the whole eigenvector matrix - just the upper half. Regarding to the link above this should be the right syntax etc.

Am I missing something? I would be very grateful for a hint.

Solution

Yes, there is a bug in LAPACKE's lapacke_cheev_work.c, in version 3.9.0, at line:

     /* Transpose output matrices */
     LAPACKE_che_trans( LAPACK_COL_MAJOR, uplo, n, a_t, lda_t, a, lda );

Indeed, LAPACKE_che_trans() is designed to transpose Hermitian matrices and makes use of their upper or lower half depending on uplo. Nevertheless, the output of cheev('V', ... ) is a matrix made by orthonormal eigenvectors of the input matrix.

The problem also occurs to other eigenvalues-related routines, such as lapacke_zheev_work.c, lapacke_zheevd_work.c

The routines lapacke_zheevr_work.c and lapacke_zheevx_work.c correctly make use of LAPACKE_zge_trans():

if( LAPACKE_lsame( jobz, 'v' ) ) {
    LAPACKE_zge_trans( LAPACK_COL_MAJOR, n, ncols_z, z_t, ldz_t, z,ldz );
}

The example provided at Intel software development tools is also likely plagued by this issue because it also makes use of LAPACKE_cheev( LAPACK_ROW_MAJOR, 'V',...), which does require transposition.

Looking at the source of LAPACKE in Lapack, version 3.8.0 is not affected by this issue, as it make use of LAPACKE_cge_trans()everywhere. Today, the problem only affects LAPACKE 3.9.0, Nov 21, 2019.

Here is a sample code to test it, chaining LAPACKE_clacpy(LAPACK_ROW_MAJOR,'U',....) and LAPACKE_cheev(LAPACK_ROW_MAJOR, 'V', 'U',...) as you do:

#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;
}

It is compiled by using

gcc main11.c -o main11b -L/home/xxxx/softs/lapack-3.9.0/lapack-3.9.0 -llapacke -llapack -lblas -lm -lgfortran -Wall

If -L/home/xxxx/softs/lapack-3.9.0/lapack-3.9.0, is removed, thus using a version of LAPACK lower than 3.9.0, it outputs:

-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 |

If the problem occurs, it prints:

-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()

The problem can be corrected by reproducing the steps performed in previous versions of 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);

I reported the issue to the github repository of Lapack.