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

学客

最近由于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

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

学客

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

学客

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

学客

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

学客

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

学客

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

学客

这样能对特征值进行从小到大顺序排序，但是特征向量不能进项从小打大排序，而是要随着特征值的顺序变化与之相应的调整，这个该怎么做到呢？非常感谢。

学客

谢谢你，这个我知道

学客

module qsort_c_module
      implicit none
      public :: QsortC
      private :: Partition
      contains
      recursive subroutine QsortC(A,B)
      real, intent(in out), dimension(:) :: A
      real, intent(in out), dimension(:,:) :: B
      integer :: iq
      if(size(A) > 1) then
      call Partition(A, iq, B)
      call QsortC(A(:iq-1),B(:,:iq-1))
      call QsortC(A(iq:),B(:,iq:))
      endif
      end subroutine QsortC
      subroutine Partition(A, marker, B)
      real, intent(in out), dimension(:) :: A
      real, intent(in out), dimension(:,:) :: B
      integer, intent(out) :: marker
      integer :: i, j
      real :: temp
      real ,dimension(1:size(A))::temb
      real :: x      ! pivot poin
      x = A(1)
      i= 0
      j= size(A) + 1
      do
          j = j-1
      do
      if (A(j) <= x) exit
      j = j-1
      end do
      i = i+1
      do
        if (A(i) >= x) exit
        i = i+1
      end do
      if (i < j) then
        ! exchange A(i) and A(j)
        temp = A(i)
        temb = B(:,i)
        A(i) = A(j)
        B(:,i)=B(:,j)
        A(j) = temp
        B(:,j)=temb
      elseif (i == j) then
        marker = i+1
        return
      else
        marker = i
        return
      endif
      end do
      end subroutine Partition
      end module qsort_c_module
      program www_fcode_cn
      use qsort_c_module
      implicit none
      integer, parameter :: r = 3
      real, dimension(1:r) :: myarray = (/0, 50, 20/)
      real, dimension(1:r,1:r) :: mb=(/1,2,3,4,5,6,7,8,9/)
      write(2,*) "原数组:", myarray
      write(2,'*') "原数组:", mb
      call QsortC(myarray,mb)
      write(2,*) "排序后:", myarray
      write(2,'*') "排序后:", mb
      end program www_fcode_cn

用这个修改后的快速排序法，排列特征值A，和特征向量B，出现了错误：
B排序输出后，出现了这样的数字：1.000000       2.000000       3.000000       7.000000   
   8.000000       9.000000       4.000000       5.000000     -1.5823616E-26
请问是哪里出错了，非常感谢