学客 发表于 2018-3-1 19:17:18

MKL求一般非对称矩阵特征值特征向量

最近由于imsl的license过期，要使用mkl计算一般非对称矩阵的特征值与特征向量，直接使用 mkl中的DGGEV函数，算出来的结果和matlab算出来的结果不一致，请大神帮忙指点:
      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

特征向量也都不一样，请大神赐教，非常感谢。

fcode 发表于 2018-3-1 23:00:44

广义特征值的求取，需要满足一定的条件，即
Ax=λBx 有非零解。

你的算例并不能满足这一条件。因此，得不到稳定的广义特征值。

你可以换一个算例，比如
A=
0.814723686393179        0.0975404049994095        0.157613081677548        0.141886338627215        0.655740699156587
0.905791937075619        0.278498218867048        0.970592781760616        0.421761282626275        0.0357116785741896
0.126986816293506        0.546881519204984        0.957166948242946        0.915735525189067        0.849129305868777
0.913375856139019        0.957506835434298        0.485375648722841        0.792207329559554        0.933993247757551
0.632359246225410        0.964888535199277        0.800280468888800        0.959492426392903        0.678735154857774

B=
0.757740130578333        0.706046088019609        0.823457828327293        0.438744359656398        0.489764395788231
0.743132468124916        0.0318328463774207        0.694828622975817        0.381558457093008        0.445586200710900
0.392227019534168        0.276922984960890        0.317099480060861        0.765516788149002        0.646313010111265
0.655477890177557        0.0461713906311539        0.950222048838355        0.795199901137063        0.709364830858073
0.171186687811562        0.0971317812358475        0.0344460805029088        0.186872604554379        0.754686681982361


这是我在 matlab 中用 rand(5) 得到的。两者都可以得到正确的，一致的特征值
Generalized eigenvalues :
(-2.49319496058144      ,0.000000000000000E+000 )
(-1.24753056286154      ,0.000000000000000E+000 )
( -6.653289900617956E-002 ,0.000000000000000E+000 )
(   1.70329737840797      ,0.000000000000000E+000 )
(   1.09592435636823      ,0.000000000000000E+000 )

学客 发表于 2018-3-4 22:11:20

非常感谢您的指点，原来是这个问题，困扰了我很长时间，我刚刚按照您的数据算了一次，特征值确实一样，但是特征向量好像不太一致，能请教您，怎么能得到一致的特征向量呢？非常感谢。

fcode 发表于 2018-3-4 22:49:55

我得到的是实数特征向量
Real*8,ALLOCATABLE::VL(:,:),VR(:,:)

建议你用 lapack95 （即 MKL 中的 95 接口）有两个好处：
1. 参数简练
2. 编译器可以帮助你检查参数类型是否错误。

use lapack95
call GGEV( RA, RB, ALPHAR, ALPHAI, BETA, VL, VR)

Eigenvectors :
-0.81514 -0.188570.210301.00000 -0.99347
0.12299 -0.91484 -0.56855 -0.146680.04254
1.000000.32993 -0.464330.25260 -0.53968
-0.57390 -0.624981.00000 -0.943501.00000
0.252621.00000 -0.25022 -0.31995 -0.41905

matlab代码是
=eig(A,B);

两者得到的一致。

学客 发表于 2018-3-5 14:58:18

好的，我用real*8 也得到与matlab相同的特征向量，但是我刚使用您说的lapack95，如果用GGEV( RA, RB, ALPHAR, ALPHAI, BETA, VL, VR),会显示编译错误。非常感谢您呢。

学客 发表于 2018-3-5 16:57:04

现在还遇到一个问题是，我用IMSL（DGVCRG）算出来的特征值和特征向量，结果都会是从最小开始排列，但是用mkl的GGEV 算出来，特征值不是按从小到大的顺序排列，这个能怎么解决呢？谢谢您。

学客 发表于 2018-3-5 20:59:17

希望得到的特征值是从小到大排列的，而且特征向量与特征值对应起来，非常感谢您

pasuka 发表于 2018-3-6 09:32:32

学客 发表于 2018-3-5 16:57
现在还遇到一个问题是，我用IMSL（DGVCRG）算出来的特征值和特征向量，结果都会是从最小开始排列，但是用mk ...
按照特征值的模，按照升序或者降序进行排序呗
本站就有很多现成代码
堆排序
http://fcode.cn/code_prof-2-1.html
快速排序
http://fcode.cn/code_prof-38-1.html

学客 发表于 2018-3-6 10:10:21

好的，谢谢呢，我去学一学

学客 发表于 2018-3-6 15:08:44

这样能对特征值进行从小到大顺序排序，但是特征向量不能进项从小打大排序，而是要随着特征值的顺序变化与之相应的调整，这个该怎么做到呢？非常感谢。
页: [1] 2
查看完整版本: MKL求一般非对称矩阵特征值特征向量