I have a problem using Lapack (DGEEV). I work with Windows XP and I use Compaq visual fortran compiler 6.6..
I wrote a numerical code for solving hyperbolic PDE that repeatedly calls the DGEEV.f subroutine to compute
eigenvalues and eigenvectors. I have the same problem using:
1)CXML libraries that are installed together with Compaq visual fortran compiler and contain both Lapack and Blas
2)lapack 3.1 and blas source codes downloaded from this site and compiled together with my codes (instead of using CXML library)
1)USING CXML libraries: I used LAPACK for a while for other problems and they always gave me correct results, but now I have a problem. While calling DGEEV
at some point of my iterations, I obtain a incorrect value of the matrix of the right eigenvectors of the following 5x5 matrix:
AI12pathMED(1,1) = 0.000000000000000E+000
AI12pathMED(2,1) = 30.3099698727331
AI12pathMED(3,1) = -30.3099698727331
AI12pathMED(4,1) = -1.201809167270153E-019
AI12pathMED(5,1) = 0.000000000000000E+000
AI12pathMED(1,2) = 0.707106781186548
AI12pathMED(2,2) = 2.612325879881627E-018
AI12pathMED(3,2) = 2.087075453413198E-019
AI12pathMED(4,2) = 7.071067811865477E-002
AI12pathMED(5,2) = 0.000000000000000E+000
AI12pathMED(1,3) = -0.707106781186548
AI12pathMED(2,3) = -1.410516712611473E-018
AI12pathMED(3,3) = 9.931016219288338E-019
AI12pathMED(4,3) = -7.071067811865477E-002
AI12pathMED(5,3) = 0.000000000000000E+000
AI12pathMED(1,4) = 0.000000000000000E+000
AI12pathMED(2,4) = 0.000000000000000E+000
AI12pathMED(3,4) = 0.000000000000000E+000
AI12pathMED(4,4) = 1.201809167270154E-018
AI12pathMED(5,4) = 0.000000000000000E+000
AI12pathMED(1,5) = 0.000000000000000E+000
AI12pathMED(2,5) = 3.195599288240487E-036
AI12pathMED(3,5) = 4.729235658520158E-037
AI12pathMED(4,5) = 1.201809167270153E-019
and DGEEV provides the following result for the matrix of the right eigenvectors:
EIGVR(1,1) = 0.000000000000000E+000
EIGVR(2,1) = 0.000000000000000E+000
EIGVR(3,1) = 0.000000000000000E+000
EIGVR(4,1) = 1.00000000000000
EIGVR(5,1) = 0.000000000000000E+000
EIGVR(1,2) = -0.150970666119249
EIGVR(2,2) = 0.698920573653124
EIGVR(3,2) = -0.698920573653125
EIGVR(4,2) = -1.509706661192487E-002
EIGVR(5,2) = 0.000000000000000E+000
EIGVR(1,3) = -0.150970666119249
EIGVR(2,3) = -0.698920573653125
EIGVR(3,3) = 0.698920573653125
EIGVR(4,3) = -1.509706661192485E-002
EIGVR(5,3) = 0.000000000000000E+000
EIGVR(1,4) = 4.933306059789007E-019
EIGVR(2,4) = -2.792685195798503E-002
EIGVR(3,4) = -2.792685195798502E-002
EIGVR(4,4) = 0.999219786573221
EIGVR(5,4) = 0.000000000000000E+000
EIGVR(1,5) = 4.933306059789008E-019
EIGVR(2,5) = -2.792685195798503E-002
EIGVR(3,5) = -2.792685195798502E-002
EIGVR(4,5) = 0.999219786573221
EIGVR(5,5) = 7.628426387804763E-258
I compared it with the results provided by the Matlab function "inv". They reads:
EIGVRm(1,1) = 0
EIGVRm(2,1) = 0
EIGVRm(3,1) = 0
EIGVRm(4,1) = 1.00000000000000
EIGVRm(5,1) = 0
EIGVRm(1,2) = -0.15097066611925
EIGVRm(2,2) = 0.69892057365312
EIGVRm(3,2) = -0.69892057365312
EIGVRm(4,2) = -0.01509706661192
EIGVRm(5,2) = 0
EIGVRm(1,3) = -0.00000000000000
EIGVRm(2,3) = -0.70708415696699
EIGVRm(3,3) = -0.70708415696699
EIGVRm(4,3) = -0.00799937076012
EIGVRm(5,3) = 0
EIGVRm(1,4) = 0.15097066611925
EIGVRm(2,4) = 0.69892057365312
EIGVRm(3,4) = -0.69892057365312
EIGVRm(4,4) = 0.01509706661192
EIGVRm(5,4) = 0
EIGVRm(1,5) = -0.00000000000000
EIGVRm(2,5) = 0.00000000000000
EIGVRm(3,5) = 0.00000000000000
EIGVRm(4,5) = -0.09950371902100
EIGVRm(5,5) = 0.99503719020999
The results are different!!!!! And then in my code I need to invert the matrix EIGVR, and I obtain a matrix with some almost infinity entries (about 1.E239). The matlab inverse matrix of EIGVRm is obviously fine since EIGVRm is ok.
2) I tried also to download the lapack 3.1 and blas source codes from this site and to compile them together with my codes instead of using CXML library.
It does not even arrive at the point were I had the problem before but it compute the wrong matrix of the eigenvectors for the following matrix:
AI12pathMED(1,1) = 0.000000000000000E+000
AI12pathMED(2,1) = -1.35233688747612
AI12pathMED(3,1) = 9.22680655780002
AI12pathMED(4,1) = 1.166596661481788E-178
AI12pathMED(5,1) = 0.000000000000000E+000
AI12pathMED(1,2) = -0.145016761150451
AI12pathMED(2,2) = -5.855729457732739E-017
AI12pathMED(3,2) = 3.708105401410622E-019
AI12pathMED(4,2) = -2.092921168058934E-163
AI12pathMED(5,2) = 0.000000000000000E+000
AI12pathMED(1,3) = 0.989429198571294
AI12pathMED(2,3) = 1.911332588533874E-016
AI12pathMED(3,3) = -3.307363042691886E-017
AI12pathMED(4,3) = 1.427971013527914E-162
AI12pathMED(5,3) = 0.000000000000000E+000
AI12pathMED(1,4) = 0.000000000000000E+000
AI12pathMED(2,4) = 0.000000000000000E+000
AI12pathMED(3,4) = 0.000000000000000E+000
AI12pathMED(4,4) = -3.054364166808208E-017
AI12pathMED(5,4) = 0.000000000000000E+000
AI12pathMED(1,5) = 0.000000000000000E+000
AI12pathMED(2,5) = -8.944365080277119E-033
AI12pathMED(3,5) = 7.540177625122472E-033
AI12pathMED(4,5) = -1.166596661481788E-178
AI12pathMED(5,5) = 0.000000000000000E+000
The matrix of the eigenvectors computed by LAPACK is again different from the one computed by matlab.Does anyone know:
1)why Lapack has problems to compute eigenvectors and Matlab not?it is strange since Matlab uses lapack, doesn't it?
2)why Lapack version 3.1 stopped before then Lapack in the CXML libraries?
3)how can I solve my problem?i.e. find a version that works?
This is the code I wrote to call the DGEEV subroutine. It computes always correct results except for the case I showed you above.
- Code: Select all
!
SUBROUTINE eigVALUEvectR(A,EIGreal,EIGimm,EIG)
!
IMPLICIT NONE
INCLUDE 'PRICE2D.DIM' ! Here I have all my common variables
!
CHARACTER JOBVL, JOBVR
INTEGER N, LDA, LDVL, LDVR, LWORK, INFO, J ,K
REAL*8 EIG(MAXVAR,MAXVAR),A(MAXVAR,MAXVAR),mat(nVAR,nVAR),
& VR(nVAR,nVAR), VL(nVAR,nVAR),WR(nVAR),WI(nVAR),WK(nVAR),
& EIGreal(MAXVAR),EIGimm(MAXVAR)
REAL*8 WORK(1:5*nVAR)
!
! Define parameter for DGEEV (eigenvector and eigenvalue) call
JOBVL = 'N'
JOBVR = 'V'
N = nVAR
LDA = nVAR
LDVL = 1
LDVR = nVAR
LWORK = 5*nvar
DO K =1,nVAR
DO J=1,nVAR
mat(K,J) = A(K,J)
ENDDO
ENDDO
!
CALL DGEEV( JOBVL, JOBVR, N, mat, LDA, WR, WI, VL, LDVL, VR,
& LDVR, WORK, LWORK, INFO )
IF (INFO.NE.0) GOTO 999
DO J = 1, N
EIGreal(j) = WR(j)
EIGimm(j) = WI(j)
ENDDO
!
DO J = 1, N
DO K = 1, LDVR
EIG(K,j) = VR(K,j)
ENDDO
ENDDO
!
RETURN
999 CONTINUE
WRITE(*,*)' IMPOSSIBLE TO FIND THE EIGENVECTORS!!'
PAUSE
STOP
!
END
thank you very much
Venexiano
