|
最近由于imsl的license过期,要使用mkl计算一般非对称矩阵的特征值与特征向量,直接使用 mkl中的DGGEV函数,算出来的结果和matlab算出来的结果不一致,请大神帮忙指点:
[Fortran] 纯文本查看 复制代码 PROGRAM MAIN
USE mkl95_PRECISION
USE mkl95_LAPACK, ONLY: GGEV
IMPLICIT NONE
CHARACTER*1 JOBVL,JOBVR
INTEGER*8 I,J,N,LDA,LDB,K,INFO,LDVL,LDVR
INTEGER*8 :: LWORK
REAL*8,ALLOCATABLE:: WORK(:)
REAL*8,ALLOCATABLE:: RA(:,:),RB(:,:)
REAL*8,ALLOCATABLE:: BETA(:),ALPHAR(:),ALPHAI(:)
COMPLEX*16,ALLOCATABLE::VL(:,:),VR(:,:)
N=5
LDA=N
LDB=N
LDVR=N
LDVL=N
LWORK=MAX(1,8*N+16)
INFO=0
ALLOCATE(RA(LDA,N),RB(LDB,N),WORK(MAX(1,LWORK)),BETA(N),
+ALPHAR(N),ALPHAI(N),VL(LDVL,N),VR(LDVR,N))
RA=0.0D0
RB=0.0D0
VR=0.0D0
VL=0.0D0
ALPHAR=0.0D0
ALPHAI=0.0D0
OPEN(5,FILE='A.INP',STATUS='OLD', ACTION='READ')! INPUT FILE
OPEN(6,FILE='A.DAT', STATUS='NEW', ACTION='WRITE') ! OUTPUT FILE
WRITE (6,*) 'GEEVX Example Program Results'
DO I=1,N
READ(5,*)(RA(I,J),J=1,N)
ENDDO
DO I=1,N
READ(5,*)(RB(I,J),J=1,N)
ENDDO
WRITE(6,*)'Matrix RA:'
DO I=1,N
WRITE(6,"(5(E13.6,1X))")RA(I,:)
ENDDO
WRITE(6,*)'Matrix RB:'
DO I=1,N
WRITE(6,"(5(E13.6,1X))")RB(I,:)
ENDDO
JOBVL='N'
JOBVR='V'
CALL DGGEV(JOBVL,JOBVR,N,RA,LDA,RB,LDB,ALPHAR,ALPHAI,BETA,VL,
+LDVL,VR,LDVR,WORK,LWORK,INFO)
WRITE(6,*)' Generalized eigenvalues : '
DO I=1,N
WRITE(6,*) '(',ALPHAR(I)/BETA(I),',',ALPHAI(I)/BETA(I),')'
ENDDO
WRITE(6,*)' Eigenvectors : '
DO I=1,N
WRITE(6,"(5(F8.5,1X,',',F8.5))") REAL(VR(I,:)),AIMAG(VR(I,:))
ENDDO
END PROGRAM MAIN
这是Fortran程序,用了5X5,A,B矩阵测试,
A =
1 2 3 4 5
6 7 8 9 10
11 12 13 14 15
16 17 18 19 20
21 22 23 24 25
B =
25 24 23 22 21
20 19 18 17 16
15 14 13 12 11
10 9 8 7 6
5 4 3 2 1
FORTRAN 算出来的特征值
Generalized eigenvalues :
( -1.00000011564385 , 0.000000000000000E+000 )
( -0.999999884356162 , 0.000000000000000E+000 )
( -3.55208454819802 , 0.000000000000000E+000 )
( -0.416408282640064 , 0.000000000000000E+000 )
( 0.284047300829184 , 0.000000000000000E+000 )
matlab的特征值:
D =
-1.0000 0 0 0 0
0 -1.0000 0 0 0
0 0 7.0565 0 0
0 0 0 1.6132 0
0 0 0 0 0.0882
特征向量也都不一样,请大神赐教,非常感谢。
|
|