dgeev_: Incorrect Eigenvalues.

[C.G.LIM]

Hi All,

Using dgeev_, I am getting correct eigenvalues for a 400x400 transposed matrix
but the eigenvalues for the original matrix are incorrect. The problem persists
even with LWORK increased to 100000 times of the required amount (i.e., 4*N).

I am using LAPACK 3.1.1 and I am calling dgeev_ from C. The eigenvalues that
I obtained for the original matrix were:

E1 = 3.64316982532374e-17 i0.00000000000000e+00
E2 = 3.64064052901039e-17 i0.00000000000000e+00
E3 = 3.63642786825317e-17 i0.00000000000000e+00

and the eigenvalues for the transposed matrix were:

E1 = 3.64316982532374e-17 i0.00000000000000e+00
E2 = 3.64064052901039e-17 i0.00000000000000e+00
E3 = 3.63642786825317e-17 i0.00000000000000e+00

The correct eigenvalues for the original 400x400 matrix are:

E1 = 4.93149e-21 + i0
E2 = 1.79213e-20 + i0
E3 = 2.86855e-20 + i0

May I know what could have gone wrong, please? Thanks!

My code is as follows:

#include<stdio.h>
#include<malloc.h>
#include<math.h>
#include"f2c.h"
#include"clapack.h"

int main(void);

int main(void)
{
double e=1.602e-19;

long int N = 400;

double *A=NULL, *B=NULL, *C=NULL, *D=NULL;
double *WR=NULL, *WI=NULL;
double *VL=NULL, *VR=NULL;
double *WORK=NULL;

char JOBVL='N', JOBVR='N';
long int LDVL=N;
long int LDVR=N;
long int LWORK=4*N;
long int INFO;

int i, row, col;

FILE *A_file, *B_file, *C_file;

double C_RowMajor[N][N], C_Transpose[N][N];

A = (double *)malloc(N*sizeof(double));
B = (double *)malloc(N*sizeof(double));
C = (double *)malloc(N*N*sizeof(double));
D = (double *)malloc(N*N*sizeof(double));

WR = (double *)malloc(N*sizeof(double));
WI = (double *)malloc(N*sizeof(double));

VL = (double *)malloc(LDVL*N*sizeof(double));
VR = (double *)malloc(LDVR*N*sizeof(double));

WORK = (double *)malloc(LWORK*sizeof(double));

for(i=0; i < N; i++)
{
if( i < 149 || i > 249 )
{
A[i] = -7.301103017227E-18;
B[i] = 1.4628975134054E-17;
}

else
{
A[i] = -9.1100330185102E-18;
B[i] = 1.822006603702E-17;
}
}

for(i=0; i < N; i++)
C[i*N+i] = B[i];

for(i=1; i < N; i++)
C[(i*N+i)-1] = A[i-1];

for(i=0; i < (N-1); i++)
C[(i*N+i)+1] = A[i+1];

C_file = fopen("C400_400.dat", "w+");

for(i=0, row=0; row < N; row++)
{
for(col=0; col < (N-1); col++)
fprintf(C_file, "%0.14e\t", C[i++]);
fprintf(C_file, "%0.14e\n", C[i++]);
}

fclose(C_file);

for(i=0, row=0; row < N; row++)
for(col=0; col < N; col++)
C_RowMajor[row][col] = C[i++];


/****** Constructing column-major array - STARTS ******/

for(i=0, col=0; col < N; col++)
for(row=0; row < N; row++)
D[i++] = C_RowMajor[row][col];

/****** Constructing column-major array - ENDS ******/


dgeev_(&JOBVL, &JOBVR, &N, D, &N, WR, WI, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);

printf("\n----------------- Original Matrix ------------------\n\nINFO = %d\n\nE1 = %0.14e i%0.14e\nE2 = %0.14e i%0.14e\nE3 = %0.14e i%0.14e\n\n----------------------------------------------------\n", INFO, WR[0], WI[0], WR[1], WI[1], WR[2], WI[2]);


/****** Transpose of Original Matrix - STARTS ******/

for(i=0, row=0; row < N; row++)
for(col=0; col < N; col++)
C_Transpose[col][row] = C[i++];

C_file = fopen("C400_400_Transpose.dat", "w+");

for(row=0; row < N; row++)
{
for(col=0; col < (N-1); col++)
fprintf(C_file, "%0.14e\t", C_Transpose[row][col]);

fprintf(C_file, "%0.14e\n", C_Transpose[row][col]);
}

fclose(C_file);

/****** Transpose of Original Matrix - ENDS ******/


/****** Constructing column-major array - STARTS ******/

for(i=0, col=0; col < N; col++)
for(row=0; row < N; row++)
D[i++] = C_Transpose[col][row];

/****** Constructing column-major array - ENDS ******/

dgeev_(&JOBVL, &JOBVR, &N, D, &N, WR, WI, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);

printf("\n--------- Transpose of the Original Matrix ---------\n\nINFO = %d\n\nE1 = %0.14e i%0.14e\nE2 = %0.14e i%0.14e\nE3 = %0.14e i%0.14e\n\n----------------------------------------------------\n\n", INFO, WR[0], WI[0], WR[1], WI[1], WR[2], WI[2]);

free(A);
free(B);
free(C);
free(D);

free(WR);
free(WI);

free(VL);
free(VR);

free(WORK);

return(0);
}

Thanks,
Cheng

[buttari]

Cheng,
I'm not sure what your problem is but I just loaded the generated files in matlab and made a quick check. Here's what I've got:

Code: Select all

>> A=load('C400_400.dat');
>> At=load('C400_400_Transpose.dat');
>> norm(At-A)

ans =

   1.8089e-18

>> norm(A'-A)

ans =

   1.8089e-18


So, A is symmetric and the eigenvalues (of both A and A') are:

Code: Select all

>> e=eig(A); e(1:3)

ans =

   1.0e-16 *

    0.3641
    0.3643
    0.3636


which is exactly what you had from the lapack routine. You should have some mistake in generating the matrix.

alfredo

[buttari]

Cheng,
I'm not sure what your problem is but I just loaded the generated files in matlab and made a quick check. Here's what I've got:

Code: Select all

>> A=load('C400_400.dat');
>> At=load('C400_400_Transpose.dat');
>> norm(At-A)

ans =

   1.8089e-18

>> norm(A'-A)

ans =

   1.8089e-18


So, A is symmetric and the eigenvalues (of both A and A') are:

Code: Select all

>> e=eig(A); e(1:3)

ans =

   1.0e-16 *

    0.3641
    0.3643
    0.3636


which is exactly what you had from the lapack routine. You should have some mistake in generating the matrix.

alfredo

[C.G.LIM]

Hi Alfredo,

Thank you very much for helping to check using Matlab. I am afraid that the matrix is not symmetrical
as it contains two pairs of unequal entries along the sub- and super-diagonal of the matrix. The two
pairs in the original matrix are:

C[148][149] = -9.11003301851020e-18
C[149][148] = -7.30110301722700e-18

and

C[148][149] = -7.30110301722700e-18
C[149][148] = -9.11003301851020e-18

I have loaded C400_400.dat and C400_400_Transpose.dat into Igor Pro 4.01A, and compute the eigenvalues
of these matrices using "MatrixEigenV MatrixEigenV /B=0 /R C400_400" and
"MatrixEigenV MatrixEigenV /B=0 /R C400_400_Transpose". The eigenvalues that I obtained for C400_400 are:

E1 = 4.93149e-21 + i0
E2 = 1.79213e-20 + i0
E3 = 2.86855e-20 + i0

whereas the eigenvalues for C400_400_Transpose are:

E1 = 3.64317e-17 + i0
E2 = 3.64064e-17 + i0
E3 = 3.63643e-17 + i0

I have checked the entries of the column-major array D that corresponds to C[148][149], C[149][148],
C[249][250], and C[250][249], and they are as follows:

D[59349] (corresponds to C[148][149]) = -9.11003301851020e-18
D[59748] (corresponds to C[149][148]) = -7.30110301722700e-18

D[99850] (corresponds to C[249][250]) = -7.30110301722700e-18
D[100249] (corresponds to C[250][249]) = -9.11003301851020e-18

I have tried multiplying all the entries by 1e18 but the problem persists.

Thanks,
Cheng

[buttari]

Hi Alfredo,

Thank you very much for helping to check using Matlab. I am afraid that the matrix is not symmetrical
as it contains two pairs of unequal entries along the sub- and super-diagonal of the matrix. The two
pairs in the original matrix are:

C[148][149] = -9.11003301851020e-18
C[149][148] = -7.30110301722700e-18

and

C[148][149] = -7.30110301722700e-18
C[149][148] = -9.11003301851020e-18

I have loaded C400_400.dat and C400_400_Transpose.dat into Igor Pro 4.01A, and compute the eigenvalues
of these matrices using "MatrixEigenV MatrixEigenV /B=0 /R C400_400" and
"MatrixEigenV MatrixEigenV /B=0 /R C400_400_Transpose". The eigenvalues that I obtained for C400_400 are:

E1 = 4.93149e-21 + i0
E2 = 1.79213e-20 + i0
E3 = 2.86855e-20 + i0

whereas the eigenvalues for C400_400_Transpose are:

E1 = 3.64317e-17 + i0
E2 = 3.64064e-17 + i0
E3 = 3.63643e-17 + i0

I have checked the entries of the column-major array D that corresponds to C[148][149], C[149][148],
C[249][250], and C[250][249], and they are as follows:

D[59349] (corresponds to C[148][149]) = -9.11003301851020e-18
D[59748] (corresponds to C[149][148]) = -7.30110301722700e-18

D[99850] (corresponds to C[249][250]) = -7.30110301722700e-18
D[100249] (corresponds to C[250][249]) = -9.11003301851020e-18

I have tried multiplying all the entries by 1e18 but the problem persists.

Thanks,
Cheng

Cheng,
I wrote the eigenvalues coming from dgeev into a file and loaded it into matlab:

Code: Select all

>> E=load('E400.dat');
>> Es=sort(E);
>> Es(1:3)

ans =

      4.89147886024384e-21
      1.78832305070559e-20
      2.86542026682222e-20


Probably you problem is just a matter of sorting. Also the check

Code: Select all

>> A=load('C400_400.dat');
>> es=sort(eig(A));
>> norm(Es-es)/norm(es)

ans =

      4.36232701005135e-15



says that the eigenvalues computed by dgeev are coorect.
hope this helps.

alfredo

[C.G.LIM]

Hi Alfredo,

Brilliant! Thank you very much for pointing out that the eigenvalues returned by dgeev_ in WR are not sorted in any order.

Cheers,
Cheng