Неправильные собственные значения с ZGEEV в LAPACK при повторении

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

В этом коде я использовал диссипативный или асимметричный c параметр lambda, который Я фиксирую как l=2.d0 и запускаю do-l oop три раза. Я должен был получить одинаковые собственные значения в каждой итерации. К моему удивлению, первые два цикла дают одинаковые наборы собственных значений, но их порядок отличается, а третий l oop генерирует совершенно другой набор собственных значений!

В чем может быть возможная ошибка?

Вывод:

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  

Примечание: есть более ранняя публикация с аналогичным возражением, но я не думайте, что я допустил ошибку, как неправильное определение RWORK, как это было сделано в этом посте. Поэтому, пожалуйста, не отмечайте мой вопрос как дубликат немедленно.

Answer 1

Проблема в том, что, цитируя страницу Згеева в

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$