Fortran Coder

标题: MKL求一般非对称矩阵特征值特征向量 [打印本页]

作者: 学客 时间: 2018-3-1 19:17
标题: 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

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

作者: fcode 时间: 2018-3-1 23:00
广义特征值的求取，需要满足一定的条件，即
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
非常感谢您的指点，原来是这个问题，困扰了我很长时间，我刚刚按照您的数据算了一次，特征值确实一样，但是特征向量好像不太一致，能请教您，怎么能得到一致的特征向量呢？非常感谢。

作者: fcode 时间: 2018-3-4 22:49
我得到的是实数特征向量
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.18857  0.21030  1.00000 -0.99347
0.12299 -0.91484 -0.56855 -0.14668  0.04254
1.00000  0.32993 -0.46433  0.25260 -0.53968
-0.57390 -0.62498  1.00000 -0.94350  1.00000
0.25262  1.00000 -0.25022 -0.31995 -0.41905

matlab代码是
[V,D]=eig(A,B);

两者得到的一致。

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

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

作者: 学客 时间: 2018-3-5 20:59
希望得到的特征值是从小到大排列的，而且特征向量与特征值对应起来，非常感谢您

作者: pasuka 时间: 2018-3-6 09: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
好的，谢谢呢，我去学一学

作者: 学客 时间: 2018-3-6 15:08
这样能对特征值进行从小到大顺序排序，但是特征向量不能进项从小打大排序，而是要随着特征值的顺序变化与之相应的调整，这个该怎么做到呢？非常感谢。

作者: pasuka 时间: 2018-3-6 15:46

学客发表于 2018-3-6 15:08
这样能对特征值进行从小到大顺序排序，但是特征向量不能进项从小打大排序，而是要随着特征值的顺序变化与之 ...

Orz。。。
特征值按照模的大小排序的同时，将对应的特征向量交换一下列或行的顺序呗

作者: 学客 时间: 2018-3-6 16:14
谢谢你，这个我知道

作者: 学客 时间: 2018-3-6 23:47
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
请问是哪里出错了，非常感谢

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)

Powered by Discuz! X3.2