Incorrect eigenvalues with ZGEEV in LAPACK when repeated

Question

I have the following code which finds the eigenvalues of a non-Hermitian complex matrix 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

In this code, I used a dissipative or asymmetric parameter lambda which I fix as l=2.d0 and run a do-loop three times. I was supposed to get the same eigenvalues in each iteration. To my surprise, the first two loops give the same sets of eigenvalues but their orderings differ, and the third loop generates a completely different set of eigenvalues!

What could be the possible error?

The output:

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

N.B.: There's an earlier post with a similar objection, but I don't think I have made a mistake like wrongly defining RWORK as it was done in that post. So please don't mark my question as a duplicate immediately.

score 4 · Accepted Answer · answered Jan 23 '20 at 21:43

The problem is that, quoting the Zgeev page at

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

A is COMPLEX*16 array, dimension (LDA,N) On entry, the N-by-N matrix A. On exit, A has been overwritten.

Emphasis is mine. Thus you need to reinitialise the whole of your matrix each time round, not just the non-zero bits. Fixing this, making your code standard conformant, and moving it into a style more consistent with modern practice I get

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$

Thanks @Ian for pointing out the flaw. Do you know when `A` is overwritten what `A` becomes on exit? The eigenvector matrix? Also, your code in the contemporary style gives `Warning: Unused dummy argument '_formal_0' at (1) [-Wunused-dummy-argument]`. — hbaromega, Jan 24 '20 at 19:32
WHo knows what it is overwritten with, and why care? The spec says it is so. That is all you really need care about. As for the warning I suspect that is an old version of gfortran not being very clever. What version are you using? I don't get that. — Ian Bush, Jan 24 '20 at 21:15
Lol! I just asked. :) I'm using Mac's Homebrew-installed GCC package: `GNU Fortran (Homebrew GCC 9.2.0_2) 9.2.0` — hbaromega, Feb 04 '20 at 11:40

Incorrect eigenvalues with ZGEEV in LAPACK when repeated

1 Answers1