重复时，LAPACK中ZGEEV的特征值不正确

Question

我有以下代码可以找到非埃尔米特复数矩阵H的特征值。

Program Eigen_nonHermitian
implicit none
integer:: Nstate,i
double precision:: t,eps,tm,tp,U
double complex:: H(6,6),E(6)
double precision:: l

! Parameters
t=1.0; eps=0.5d0; U = 2.d0
Nstate=6

! Initialize H
H=dcmplx(0.d0,0.d0)

l=2.d0

do i=1,3
    ! Construct H-matrix (size: 6 x 6)
    tm=t-l; tp=t+l

    !H(1,1)=2.d0*eps; H(3,3)=H(1,1); H(4,4)=H(1,1); H(6,6)=H(1,1)
    !H(2,2)=2.d0*eps+U; H(5,5)=H(2,2)
    !H(2,3)=tm; H(3,5)=tm
    !H(3,2)=tp; H(5,3)=tp
    !H(2,4)=-tm; H(4,5)=-tm
    !H(4,2)=-tp; H(5,4)=-tp

    H(1,1)=dcmplx(2.d0*eps,0.d0); H(3,3)=H(1,1); H(4,4)=H(1,1); H(6,6)=H(1,1)
    H(2,2)=dcmplx(2.d0*eps+U,0.d0); H(5,5)=H(2,2)
    H(2,3)=dcmplx(tm,0.d0); H(3,5)=dcmplx(tm,0.d0)
    H(3,2)=dcmplx(tp,0.d0); H(5,3)=dcmplx(tp,0.d0)
    H(2,4)=dcmplx(-tm,0.d0); H(4,5)=dcmplx(-tm,0.d0)
    H(4,2)=dcmplx(-tp,0.d0); H(5,4)=dcmplx(-tp,0.d0)

    ! Diagonalize H 
    call CEigen(H,E,Nstate)


    print*,'<== Dissipative parameter, lambda=',l

end do


End Program Eigen_nonHermitian

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Subroutine CEigen(A,ev,N)
        implicit none
        integer:: i
    integer:: N,LDA,LDVL,LDVR,LWORK,INFO
    double complex:: A(N,N),ev(N)
    double complex:: vecL(N,N), vecR(N,N), WORK(2*N)
    double precision:: RWORK(2*N)

    LDA=N; LDVL=N; LDVR=N

        LWORK=2*N 
        RWORK=2*N

    call ZGEEV('V', 'V', N, A, LDA, ev, vecL, LDVL, vecR, LDVR, &
         & WORK, LWORK, RWORK, info)
    ! Syntax:  ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
                 !     WORK, LWORK, RWORK, INFO )

    print*,'Eigenvalues are:' 
    do i=1, N
           write(*,'(f8.3,a,f8.3,a)') dreal(ev(i)),' +', & 
                & dimag(ev(i)),'i'
    enddo
        print*,'info=',info

end Subroutine CEigen

在此代码中，我使用了耗散性或非对称性参数lambda，将其固定为l=2.d0，并运行了3次do-loop。我应该在每次迭代中得到相同的特征值。令我惊讶的是，前两个循环给出了相同的特征值集，但是它们的顺序不同，而第三个循环产生了完全不同的特征值集！

可能是什么错误？

输出：

Eigenvalues are:
   2.000 +   3.317i
   3.000 +  -0.000i
   2.000 +  -3.317i
   1.000 +   0.000i
   1.000 +   0.000i
   1.000 +   0.000i
 info=           0
 <== Dissipative parameter, lambda=   2.0000000000000000     
 Eigenvalues are:
   1.000 +   0.000i
   2.000 +  -3.317i
   3.000 +   0.000i
   2.000 +   3.317i
   1.000 +   0.000i
   1.000 +   0.000i
 info=           0
 <== Dissipative parameter, lambda=   2.0000000000000000     
 Eigenvalues are:
   1.000 +   0.000i
   2.140 +   3.421i
   2.180 +  -3.329i
   3.000 +   0.000i
   0.680 +  -0.092i
   1.000 +   0.000i
 info=           0
 <== Dissipative parameter, lambda=   2.0000000000000000  

NB .:有一个earlier post也有类似的异议，但我不认为我犯了像错误定义RWORK一样的错误。因此，请不要立即将我的问题标记为重复。

Answer 1

问题是，引用Zgeev页面位于

http://www.netlib.org/lapack/explore-html/db/d55/group__complex16_g_eeigen_ga0eb4e3d75621a1ce1685064db1ac58f0.html

A是COMPLEX * 16数组，尺寸（LDA，N）进入时，N×N矩阵A。在出口处，A已被覆盖。

强调是我的。因此，您需要在每个时间重新初始化整个矩阵，而不仅仅是非零位。解决此问题，使您的代码符合标准，并将其转换为更符合现代实践的样式，我得到

ian@eris:~/work/stack$ cat eig.f90
Program Eigen_nonHermitian
  implicit none
  Integer, Parameter :: wp = Selected_real_kind( 12, 70 )
  integer:: Nstate,i
  Real( wp ) :: t,eps,tm,tp,U
  Complex( wp ) :: H(6,6),E(6)
  Real( wp ) :: l

  ! Parameters
  t=1.0_wp; eps=0.5_wp; U = 2.0_wp
  Nstate=6

  ! Initialize H

  l=2.0_wp

  do i=1,3
     H=Cmplx(0.0_wp,0.0_wp, Kind = Kind( H ) )
     ! Construct H-matrix (size: 6 x 6)
     tm=t-l; tp=t+l

     !H(1,1)=2.d0*eps; H(3,3)=H(1,1); H(4,4)=H(1,1); H(6,6)=H(1,1)
     !H(2,2)=2.d0*eps+U; H(5,5)=H(2,2)
     !H(2,3)=tm; H(3,5)=tm
     !H(3,2)=tp; H(5,3)=tp
     !H(2,4)=-tm; H(4,5)=-tm
     !H(4,2)=-tp; H(5,4)=-tp

     H(1,1)=Cmplx(2._wp*eps,0._wp, Kind = Kind( H ) ); H(3,3)=H(1,1); H(4,4)=H(1,1); H(6,6)=H(1,1)
     H(2,2)=Cmplx(2._wp*eps+U,0._wp, Kind = Kind( H ) ); H(5,5)=H(2,2)
     H(2,3)=Cmplx(tm,0._wp, Kind = Kind( H ) ); H(3,5)=Cmplx(tm,0._wp, Kind = Kind( H ) )
     H(3,2)=Cmplx(tp,0._wp, Kind = Kind( H ) ); H(5,3)=Cmplx(tp,0._wp, Kind = Kind( H ) )
     H(2,4)=Cmplx(-tm,0._wp, Kind = Kind( H ) ); H(4,5)=Cmplx(-tm,0._wp, Kind = Kind( H ) )
     H(4,2)=Cmplx(-tp,0._wp, Kind = Kind( H )); H(5,4)=Cmplx(-tp,0._wp, Kind = Kind( H ) )

     ! Diagonalize H 
     call CEigen(H,E,Nstate)


     print*,'<== Dissipative parameter, lambda=',l

  end do

Contains

  Subroutine CEigen(A,ev,N)
    implicit none
    integer:: i
    integer:: N,LDA,LDVL,LDVR,LWORK,INFO
    complex( wp ):: A(N,N),ev(N)
    complex( wp ):: vecL(N,N), vecR(N,N), WORK(2*N)
    Real( wp ):: RWORK(2*N)

    LDA=N; LDVL=N; LDVR=N

    LWORK=2*N 
    RWORK=2*N

    call ZGEEV('V', 'V', N, A, LDA, ev, vecL, LDVL, vecR, LDVR, &
         & WORK, LWORK, RWORK, info)
    ! Syntax:  ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
    !     WORK, LWORK, RWORK, INFO )

    print*,'Eigenvalues are:' 
    do i=1, N
       write(*,'(f8.3,a,f8.3,a)') Real(ev(i), Kind = Kind( ev ) ),' +', & 
            & Aimag(ev(i)),'i'
    enddo
    print*,'info=',info

  end Subroutine CEigen

End Program Eigen_nonHermitian

ian@eris:~/work/stack$ gfortran --version
GNU Fortran (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ian@eris:~/work/stack$ gfortran -std=f2008 -Wall -Wextra -fcheck=all eig.f90 -llapack
ian@eris:~/work/stack$ ./a.out
 Eigenvalues are:
   2.000 +   3.317i
   3.000 +   0.000i
   2.000 +  -3.317i
   1.000 +   0.000i
   1.000 +   0.000i
   1.000 +   0.000i
 info=           0
 <== Dissipative parameter, lambda=   2.0000000000000000     
 Eigenvalues are:
   2.000 +   3.317i
   3.000 +   0.000i
   2.000 +  -3.317i
   1.000 +   0.000i
   1.000 +   0.000i
   1.000 +   0.000i
 info=           0
 <== Dissipative parameter, lambda=   2.0000000000000000     
 Eigenvalues are:
   2.000 +   3.317i
   3.000 +   0.000i
   2.000 +  -3.317i
   1.000 +   0.000i
   1.000 +   0.000i
   1.000 +   0.000i
 info=           0
 <== Dissipative parameter, lambda=   2.0000000000000000     
ian@eris:~/work/stack$