Question about LAPACK in FORTRAN 90 (about ZGEEV)

[EM_waves]

I have some question about zgeev subroutine in LAPCK to trouble you.

I am writing FORTRAN 90 code. zgeev subroutine written in F77 is not compatible with my FORTRAN 90 code (zgeev subroutine written in F77 does not return the correct eigenvectors and eigenvalues. Even worse when zgeev subroutine is called the second times, I encountered error: servere(161): Program Exception – array bounds exceeded at line 167 in ZGEBAK.f). I found that if the arguments of zgeev were not allocatable matrices, I did not encounter the error. So I assume this is due to the zgeev subroutine written in F77 is not compatible with my FORTRAN 90 code.

Is the LAPACK written in FORTRAN 90 available? I found the source codes of LAPACK in FORTRAN 95. But when I tried to compile the source codes to generate library, I had more than 700 errors.

What should I do? Thank you very much for your help.

The following is my FORTRAN 90 code that has problem with ZGEEV in FORTRAN 77. When all matrices are not allocatable matices, I do not have the error. But I need to use allocatable matrices. Thanks again.

program test_eig

implicit none
INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0))
INTEGER, PARAMETER :: DP = KIND(1.0D0)
complex(DPC), dimension(:,:), allocatable :: VL,VR, A
complex(DPC), dimension(:), allocatable :: W, WORK
real(DP), dimension(:), allocatable :: RWORK
integer :: INFO, N

N=6

allocate(VL(1,1))
allocate(VR(N,N), A(N,N))
allocate(W(N), WORK(N))
allocate(RWORK(N*2))
VL=0
VR=0
W=0
WORK=0
RWORK=0

A=0
A(1,1)=(0.948275170054075,1.13160317511797)
A(2,1)=(-0.849373459407665,-1.90244480766304)
A(3,1)=(6.97177713356470,-5.50215064828530)

A(1,2)=(-16.8336381894322,29.0327833705286)
A(2,2)=(-0.481557522697477,8.350384821317058E-002)
A(3,2)=(-8.62609828223245,-32.4324577138091)

A(1,3)=(5.50215064828525,6.97177713356471)
A(2,3)=(-0.744633891424639,1.94582935723871)
A(3,3)=(0.948275170054117,1.13160317511794)

A(4,4)=(-6.69287772055732,7.050000131130217E-002)
A(5,4)=(-0.720175500523152,-1.61830643033964)
A(6,4)=0

A(4,5)=(-0.635074470866619,1.65355643099529)
A(5,5)=(-0.409499942779541,7.050000131130217E-002)
A(6,5)=(-0.720175500523152,-1.61830643033964)


A(5,6)=(-0.635074470866619,1.65355643099529)
A(6,6)=(-6.69287772055732,7.050000131130217E-002)
A(4,6)=0

call ZGEEV('N','V',N,A,N,W,VL,1,VR, N, WORK, N*2, RWORK, INFO)

A=0
A(1,1)=(-9.40337766333686,0.141000002622604)
A(2,2)=(-3.11999988555908,0.141000002622604)
A(3,3)=(-9.40337766333686,0.141000002622604)
A(4,4)=(-9.40337766333686,0.141000002622604)
A(5,5)=(-3.11999988555908,0.141000002622604)
A(6,6)=(-9.40337766333686,0.141000002622604)
call ZGEEV('N','V',N,A,N,W,VL,1,VR, N, WORK, N*2, RWORK, INFO)


end program test_eig

[sven]

Dear EM_waves,

It looks to me as though you have not allocated enough space to WORK. I think that you should have

allocate(W(N), WORK(2*N))

Note that if you are expecting to solve problems for large values of N then, for efficiency, the length of WORK should be larger; something like (1+NB)*N, where NB is the block size. Taking NB = 64 is probably reasonable.

Best wishes,

Sven Hammarling.

[EM_waves]

Thanks a lot Sven!

Now the problem is solved. But I have a further question. I compared the results with Matlab. I got the same eigenvalues in the same order. But the order of some elements of the eigenvectors I obtain from my Fortran code are different from that of Matlab. Please notice that
VR(1,1)
VR(2,1)
VR(3,1)
are different. What is wrong this time? I only changed
allocate(W(N), WORK(2*N)) in the Fortran code I posted in my last thread. Thank you very much.

The eigenvectors I obtained from Fortran are,

VR(1,1) (0.695524450618616,0.000000000000000E+000)
VR(2,1) (7.947120078562600E-002,0.161789393587226)
VR(3,1) (-0.491810056057660,0.491810052480225)
VR(4,1) (0.000000000000000E+000,0.000000000000000E+000)
VR(5,1) (0.000000000000000E+000,0.000000000000000E+000)
VR(6,1) (0.000000000000000E+000,0.000000000000000E+000)
VR(1,2) (0.499999992503837,0.500000004879650)
VR(2,2) (1.792354965359064E-009,1.763453806851868E-009)
VR(3,2) (0.707106783036702,0.000000000000000E+000)
VR(4,2) (0.000000000000000E+000,0.000000000000000E+000)
VR(5,2) (0.000000000000000E+000,0.000000000000000E+000)
VR(6,2) (0.000000000000000E+000,0.000000000000000E+000)
VR(1,3) (0.669720526364891,0.000000000000000E+000)
VR(2,3) (-0.150655228044142,-0.283287561731985)
VR(3,3) (-0.473563874656229,0.473563969086572)
VR(4,3) (0.000000000000000E+000,0.000000000000000E+000)
VR(5,3) (0.000000000000000E+000,0.000000000000000E+000)
VR(6,3) (0.000000000000000E+000,0.000000000000000E+000)
VR(1,4) (0.000000000000000E+000,0.000000000000000E+000)
VR(2,4) (0.000000000000000E+000,0.000000000000000E+000)
VR(3,4) (0.000000000000000E+000,0.000000000000000E+000)
VR(4,4) (-8.496613311564774E-002,0.217533941496527)
VR(5,4) (0.943885313644129,0.000000000000000E+000)
VR(6,4) (-9.373959597293809E-002,-0.213899844212461)
VR(1,5) (0.000000000000000E+000,0.000000000000000E+000)
VR(2,5) (0.000000000000000E+000,0.000000000000000E+000)
VR(3,5) (0.000000000000000E+000,0.000000000000000E+000)
VR(4,5) (0.707106795034171,0.000000000000000E+000)
VR(5,5) (1.214143197981906E-016,3.242507729843877E-016)
VR(6,5) (0.499999982327747,-0.499999998088757)
VR(1,6) (0.000000000000000E+000,0.000000000000000E+000)
VR(2,6) (0.000000000000000E+000,0.000000000000000E+000)
VR(3,6) (0.000000000000000E+000,0.000000000000000E+000)
VR(4,6) (0.667427719010721,0.000000000000000E+000)
VR(5,6) (0.132567810552445,0.302500066598763)

What I obtained from matlab is

VR(1,1) -0.4918 - 0.4918i
VR(2,1) 0.0582 - 0.1706i
VR(3,1) 0.6955
VR(4,1) 0
VR(5,1) 0
VR(6,1) 0
VR(1,2) 0.7071
VR(2,2) -0.0000 + 0.0000i
VR(3,2) 0.5000 - 0.5000i
VR(4,2) 0
VR(5,2) 0
VR(6,2) 0
VR(1,3) 0.6697
VR(2,3) -0.1507 - 0.2833i
VR(3,3) -0.4736 + 0.4736i
VR(4,3) 0
VR(5,3) 0
VR(6,3) 0
VR(1,4) 0
VR(2,4) 0
VR(3,4) 0
VR(4,4) -0.0850 + 0.2175i
VR(5,4) 0.9439
VR(6,4) -0.0937 - 0.2139i
VR(1,5) 0
VR(2,5) 0
VR(3,5) 0
VR(4,5) 0.7071
VR(5,5) -0.0000 + 0.0000i
VR(6,5) 0.5000 - 0.5000i
VR(1,6) 0
VR(2,6) 0
VR(3,6) 0
VR(4,6) -0.4719 - 0.4719i
VR(5,6) 0.1202 - 0.3076i
VR(6,6) 0.6674

Dear EM_waves,

It looks to me as though you have not allocated enough space to WORK. I think that you should have

allocate(W(N), WORK(2*N))

Note that if you are expecting to solve problems for large values of N then, for efficiency, the length of WORK should be larger; something like (1+NB)*N, where NB is the block size. Taking NB = 64 is probably reasonable.

Best wishes,

Sven Hammarling.

[Julien Langou]

Hello,
I think you messed up your Matlab code.
With the following code I do have the same VR has you got from your Fortran code.
Julien

Code: Select all

A=zeros(6);
A(1,1)=0.948275170054075+sqrt(-1)*1.13160317511797;
A(2,1)=-0.849373459407665+sqrt(-1)*-1.90244480766304;
A(3,1)=6.97177713356470+sqrt(-1)*-5.50215064828530;
A(1,2)=-16.8336381894322+sqrt(-1)*29.0327833705286;
A(2,2)=-0.481557522697477+sqrt(-1)*8.350384821317058E-002;
A(3,2)=-8.62609828223245+sqrt(-1)*-32.4324577138091;
A(1,3)=5.50215064828525+sqrt(-1)*6.97177713356471;
A(2,3)=-0.744633891424639+sqrt(-1)*1.94582935723871;
A(3,3)=0.948275170054117+sqrt(-1)*1.13160317511794;
A(4,4)=-6.69287772055732+sqrt(-1)*7.050000131130217E-002;
A(5,4)=-0.720175500523152+sqrt(-1)*-1.61830643033964;
A(4,5)=-0.635074470866619+sqrt(-1)*1.65355643099529;
A(5,5)=-0.409499942779541+sqrt(-1)*7.050000131130217E-002;
A(6,5)=-0.720175500523152+sqrt(-1)*-1.61830643033964;
A(5,6)=-0.635074470866619+sqrt(-1)*1.65355643099529;
A(6,6)=-6.69287772055732+sqrt(-1)*7.050000131130217E-002;
[v,d]=eig(A);
v


[EM_waves]

Thank you so much Julien. That really helps. I found out the problem.

In matlab, the following code,
A=zeros(6);
A(1,1)=0.948275170054075+sqrt(-1)*1.13160317511797;
A(2,1)=-0.849373459407665+sqrt(-1)*-1.90244480766304;
A(3,1)=6.97177713356470+sqrt(-1)*-5.50215064828530;
A(1,2)=-16.8336381894322+sqrt(-1)*29.0327833705286;
A(2,2)=-0.481557522697477+sqrt(-1)*8.350384821317058E-002;
A(3,2)=-8.62609828223245+sqrt(-1)*-32.4324577138091;
A(1,3)=5.50215064828525+sqrt(-1)*6.97177713356471;
A(2,3)=-0.744633891424639+sqrt(-1)*1.94582935723871;
A(3,3)=0.948275170054117+sqrt(-1)*1.13160317511794;
A(4,4)=-6.69287772055732+sqrt(-1)*7.050000131130217E-002;
A(5,4)=-0.720175500523152+sqrt(-1)*-1.61830643033964;
A(4,5)=-0.635074470866619+sqrt(-1)*1.65355643099529;
A(5,5)=-0.409499942779541+sqrt(-1)*7.050000131130217E-002;
A(6,5)=-0.720175500523152+sqrt(-1)*-1.61830643033964;
A(5,6)=-0.635074470866619+sqrt(-1)*1.65355643099529;
A(6,6)=-6.69287772055732+sqrt(-1)*7.050000131130217E-002;
[v1,d1]=eig(A);

I got eigenvectors,
VR(1,1) 0.695524450254131 + 0i
VR(2,1) 0.079471200901875 + 0.16178938909766i
VR(3,1) -0.49181005525574 + 0.491810055255741i
VR(4,1) 0 + 0i
VR(5,1) 0 + 0i
VR(6,1) 0 + 0i
VR(1,2) 0.49999999999999 + 0.500000000000007i
VR(2,2) 0.000000000000002 + 0.000000000000001i
VR(3,2) 0.70710678118655 + 0i
VR(4,2) 0 + 0i
VR(5,2) 0 + 0i
VR(6,2) 0 + 0i
VR(1,3) 0.669720524244207 + 0i
VR(2,3) -0.150655217423941 - 0.283287564631902i
VR(3,3) -0.473563924192885 + 0.473563924192891i
VR(4,3) 0 + 0i
VR(5,3) 0 + 0i
VR(6,3) 0 + 0i
VR(1,4) 0 + 0i
VR(2,4) 0 + 0i
VR(3,4) 0 + 0i
VR(4,4) -0.084966130230682 + 0.217533938347771i
VR(5,4) 0.943885313351714 + 0i
VR(6,4) -0.093739596086631 - 0.213899849801219i
VR(1,5) 0 + 0i
VR(2,5) 0 + 0i
VR(3,5) 0 + 0i
VR(4,5) 0.707106781186548 + 0i
VR(5,5) -0 - 0i
VR(6,5) 0.499999999999999 - 0.500000000000001i
VR(1,6) 0 + 0i
VR(2,6) 0 + 0i
VR(3,6) 0 + 0i
VR(4,6) -0.471942656675855 - 0.471942656675857i
VR(5,6) 0.120160253714588 - 0.307639445887849i
VR(6,6) 0.667427705733387 + 0i

When I changed the matlab code to
A=zeros(6);
A(1,1)=0.948274308785464 + 1.131602862739101i;
A(2,1)=-0.849373476348601 - 1.902444850455919i;
A(3,1)=6.971776357250898 - 5.502150312733953i;
A(1,2)=-16.833635696567036 + 29.032780759784607i;
A(2,2)=-0.481557590470529 + 0.083503846874571i ;
A(3,2)=-8.626098198879543 - 32.4324541050125i;
A(1,3)=5.502150312733981 + 6.971776357250896i;
A(2,3)=-0.74463390970472 + 1.945829399476895i;
A(3,3)=0.948274308785446 + 1.131602862739129i;
A(4,4)=-6.692877777777779 + 0.0705i;
A(5,4)=-0.720175514036826 - 1.618306465203068i;
A(4,5)=-0.635074485963174 + 1.653556465203068i;
A(5,5)=-0.4095 + 0.0705i;
A(6,5)=-0.720175514036826 - 1.618306465203068i;
A(5,6)=-0.635074485963174 + 1.653556465203068i;
A(6,6)=-6.692877777777779 + 0.0705i;
[v2,d2]=eig(A);

I got eigenvectors
VR(1,1) -0.49181005444791 - 0.491810054447911i
VR(2,1) 0.05820775366156 - 0.170597006976817i
VR(3,1) 0.695524449111685 + 0i
VR(4,1) 0 + 0i
VR(5,1) 0 + 0i
VR(6,1) 0 + 0i
VR(1,2) 0.500000000000005 + 0.499999999999995i
VR(2,2) -0.000000000000001 + 0i
VR(3,2) 0.707106781186548 + 0i
VR(4,2) 0 + 0i
VR(5,2) 0 + 0i
VR(6,2) 0 + 0i
VR(1,3) 0.669720518999101 + 0i
VR(2,3) -0.150655223152972 - 0.283287586385072i
VR(3,3) -0.473563920484037 + 0.473563920484038i
VR(4,3) 0 + 0i
VR(5,3) 0 + 0i
VR(6,3) 0 + 0i
VR(1,4) 0 + 0i
VR(2,4) 0 + 0i
VR(3,4) 0 + 0i
VR(4,4) -0.084966131689905 + 0.217533941458682i
VR(5,4) 0.943885311655081 + 0i
VR(6,4) -0.09373959725455 - 0.213899853032792i
VR(1,5) 0 + 0i
VR(2,5) 0 + 0i
VR(3,5) 0 + 0i
VR(4,5) 0.707106781186548 + 0i
VR(5,5) -0 + 0i
VR(6,5) 0.499999999999999 - 0.5i
VR(1,6) 0 + 0i
VR(2,6) 0 + 0i
VR(3,6) 0 + 0i
VR(4,6) -0.47194265582754 - 0.47194265582754i
VR(5,6) 0.120160255778242 - 0.307639450287343i
VR(6,6) 0.667427704533686 + 0i

The eigenproblem is very sensitive to the matrix itself. The difference of the two matrices is at 6 digits.

[Julien Langou]

Hello,

hop hop, I did not exactly say that. Your eigenproblem is extremely well conditionned. What I have said is that if you give exactly the same input matrix then MATLAB and LAPACK returns exactly the same answer, which makes some sense since MATLAB uses LAPACK subroutines.

Now why when you modify the data but keeping up to 6 digits, do you get that big difference in your eigenvectors?

Keep in mind that an eigenvector V is defined up to a scaling coefficient.

So when the eigenvectors are normalized ( ||v||_2 = 1, normalization is done by LAPACK and MATLAB) this means that two eigenvectors output from two different eigensolvers can differ by any multiplication by e^{i.theta}, so when e^{i.theta} = +/- 1 you can check this easily with your eyes but otherwise it is not obvious to check.

In Matlab, name A1 your first matrix, A2 the second. Then:

Code: Select all

[v1,d1]=eig(A1);
[v2,d2]=eig(A2);
s=zeros(6,1);
s(1)=1;
s(2)=v1(1,2)/v2(1,2);
s(3)=v1(1,3)/v2(1,3);
s(4)=1;
s(5)=v1(4,5)/v2(4,5);
s(6)=v1(4,6)/v2(4,6);
norm(A1-A2) / norm(A1)
norm(d1-d2)/norm(d1)
norm( v1/diag(s) - v2 )


And you should get
|| A1 - A2 || / || A1 || = 1e-7
|| D1 - D2 || / || D1 || = 1e-7
|| V1 *S - V2 || / || V1 || = 2e-8
where S is a diagonal matrix with diagonal elements of absolute value 1.

-j

[EM_waves]

Thank you again Julien. I really learn something from you.

I calculated the following in matlab

A1=zeros(6);
A1(1,1)=0.948275170054075+sqrt(-1)*1.13160317511797;
A1(2,1)=-0.849373459407665+sqrt(-1)*-1.90244480766304;
A1(3,1)=6.97177713356470+sqrt(-1)*-5.50215064828530;
A1(1,2)=-16.8336381894322+sqrt(-1)*29.0327833705286;
A1(2,2)=-0.481557522697477+sqrt(-1)*8.350384821317058E-002;
A1(3,2)=-8.62609828223245+sqrt(-1)*-32.4324577138091;
A1(1,3)=5.50215064828525+sqrt(-1)*6.97177713356471;
A1(2,3)=-0.744633891424639+sqrt(-1)*1.94582935723871;
A1(3,3)=0.948275170054117+sqrt(-1)*1.13160317511794;
A1(4,4)=-6.69287772055732+sqrt(-1)*7.050000131130217E-002;
A1(5,4)=-0.720175500523152+sqrt(-1)*-1.61830643033964;
A1(4,5)=-0.635074470866619+sqrt(-1)*1.65355643099529;
A1(5,5)=-0.409499942779541+sqrt(-1)*7.050000131130217E-002;
A1(6,5)=-0.720175500523152+sqrt(-1)*-1.61830643033964;
A1(5,6)=-0.635074470866619+sqrt(-1)*1.65355643099529;
A1(6,6)=-6.69287772055732+sqrt(-1)*7.050000131130217E-002;
[v1,d1]=eig(A1);

A2=zeros(6);
A2(1,1)=0.948274308785464 + 1.131602862739101i;
A2(2,1)=-0.849373476348601 - 1.902444850455919i;
A2(3,1)=6.971776357250898 - 5.502150312733953i;
A2(1,2)=-16.833635696567036 + 29.032780759784607i;
A2(2,2)=-0.481557590470529 + 0.083503846874571i ;
A2(3,2)=-8.626098198879543 - 32.4324541050125i;
A2(1,3)=5.502150312733981 + 6.971776357250896i;
A2(2,3)=-0.74463390970472 + 1.945829399476895i;
A2(3,3)=0.948274308785446 + 1.131602862739129i;
A2(4,4)=-6.692877777777779 + 0.0705i;
A2(5,4)=-0.720175514036826 - 1.618306465203068i;
A2(4,5)=-0.635074485963174 + 1.653556465203068i;
A2(5,5)=-0.4095 + 0.0705i;
A2(6,5)=-0.720175514036826 - 1.618306465203068i;
A2(5,6)=-0.635074485963174 + 1.653556465203068i;
A2(6,6)=-6.692877777777779 + 0.0705i;
[v2,d2]=eig(A2);

ratio=abs(v1./v2);

I expect that the elements of ratio should be all 1 except the elements corresponding to zero elements of v2. I got ratio as in the following

ratio =

1.0000 1.0000 1.0000 NaN NaN NaN
1.0000 1.6047 1.0000 NaN NaN NaN
1.0000 1.0000 1.0000 NaN NaN NaN
NaN NaN NaN 1.0000 1.0000 1.0000
NaN NaN NaN 1.0000 1.1802 1.0000
NaN NaN NaN 1.0000 1.0000 1.0000

1.6047 and 1.1802 are all right since they corresponds to almost zero elements in eigenvector matrices.

[Julien Langou]

ok good you almost got it.
So if you look carefully those 1.??? 1.??? corresponds to very small entry in v1 and v2 so they are meaning less and so you can/should discard them.

So if you discard anything smaller than 1e-13 (for example):

Code: Select all

n=6;ratio=zeros(n);for j=1:n, for i=1:n,
if (abs(v1(i,j)) >1e-13 ), ratio(i,j)=v1(i,j)/v2(i,j);
end;end;end

ratio =
   1.0000             0.7071 + 0.7071i  -0.7071 + 0.7071i        0                  0                  0         
   1.0000 - 0.0000i        0            -0.7071 + 0.7071i        0                  0                  0         
   1.0000 - 0.0000i   0.7071 + 0.7071i  -0.7071 + 0.7071i        0                  0                  0         
        0                  0                  0             1.0000 - 0.0000i   0.7071 - 0.7071i  -0.7071 - 0.7071i
        0                  0                  0             1.0000                  0            -0.7071 - 0.7071i
        0                  0                  0             1.0000 - 0.0000i   0.7071 - 0.7071i  -0.7071 - 0.7071i



-j

[EM_waves]

Thanks