Is biorthogonalization of left and right eigenvectors using ZGEEV in lapack guaranteed?

Question

Let's say I want to diagonalize a complex non-Hermitian matrix of the form H = H0 + iV, where H0 and V are Hermitian matrices. Let R and L are matrices containing right and left eigenvectors respectively as the columns. Then conjugate(transpose(L)).R = I should be satisfied according to biorthogonalization condition. But L and R as obtained using ZGEEV does not seem to satisfy this condition. Please look into the following fortran90 code:

PROGRAM NonHerm

implicit none

INTEGER,Parameter :: m = 20,   avgstep = 1 

integer :: i,j,k,q,p

real*8, allocatable ::  eig(:)

complex*16, allocatable :: h(:,:), W(:), evl(:,:), evr(:,:)

complex*16 :: im, t, g, h2(m,m), norm

allocate(h(m,m))

allocate(W(m))

allocate(evl(m,m))

allocate(evr(m,m))

im = (0.0d0,1.d0);
t = 1.d0;
g = 0.0d0

call ham(m,h,t,g,im)

call compl_diag_lr(m,h,W,evr,evl)

do j = 1, m

    norm = 0.d0

    do k = 1, m

        norm = norm + conjg(evl(k,j))*evr(k,j)

    end do

    evr(:,j) = evr(:,j)/norm

end do

h2 = matmul(conjg(transpose(evl)), evr)

do i = 1, m

    print*, (h2(i,j), j = 1, m)  !h2 should be an Indentity matrix

end do

deallocate(h)

deallocate(W)

deallocate(evl)

deallocate(evr) 

END PROGRAM

subroutine ham(d,h,t,g,im)

implicit none

integer :: i,j,k

integer :: d

complex*16 :: im, t, g

complex*16 :: h(d,d)

h = 0.d0

do i = 1, d

    j = mod(i,d) + 1

    h(i,j) = t

    h(j,i) = conjg(t)

    h(i,i) = im*g

end do

end subroutine

subroutine compl_diag_lr(d,h,W,VR,VL)

implicit none

      integer :: d

      INTEGER :: INFO,LWORK

      COMPLEX*16, DIMENSION(2*d) :: WORK

      REAL*8, DIMENSION(2*d) :: RWORK

      COMPLEX*16, DIMENSION(d,d) :: VR, h

      COMPLEX*16, DIMENSION(d,d) :: VL

      COMPLEX*16, DIMENSION(d) :: W   

LWORK = 2*d

CALL ZGEEV( 'V', 'V', d, h, d, W, VL, d, VR, d, WORK,LWORK, RWORK, INFO ) 

end subroutine

Ian Bush · Accepted Answer

Your matrix has degenerate eigenvalues. Biorthogonality is not guaranteed between evecs corresponding to such evals. It is these cases that are non-zero, so there is nothing wrong with LAPACK. Note biorthogonality can be enforced in such cases, but the LAPACK documentation does not say it will do such a thing, so you can't depend upon it - as you see here. If you require it you will have to write your own code to do it.

Note also your code is not standard Fortran - complex*16 and similar is not and has never been part of Fortran and should not be used. Please learn about Fortran kinds . Also external subroutines should be deprecated, use contained or, best, module subprograms instead. Here is a standard conforming version of your program with more modern programming style, and easier to read output:

ijb@ijb-Latitude-5410:~/work/stack$ cat zgeev.f90
Module matrix_routines

Contains

  Subroutine ham(h,t,g)

    Use, Intrinsic :: iso_fortran_env, Only : wp => real64
    
    Implicit None

    Complex( wp ), Dimension( :, : ), Intent(   Out ) :: h(:,:)
    Complex( wp ),                    Intent( In    ) :: t, g

    Integer :: i,j

    Integer :: d

    d = Size( h, Dim = 1 )

    h = 0.0_wp

    Do i = 1, d

       j = Mod(i,d) + 1

       h(i,j) = t

       h(j,i) = Conjg(t)

       h(i,i) = ( 0.0_wp, 1.0_wp ) * g

    End Do

  End Subroutine ham

  Subroutine compl_diag_lr(h,W,VR,VL)

    Use, Intrinsic :: iso_fortran_env, Only : wp => real64

    Implicit None

    Complex( wp ), Dimension(:,:), Intent( InOut ) :: h
    Complex( wp ), Dimension(:  ), Intent(   Out ) :: W
    Complex( wp ), Dimension(:,:), Intent(   Out ) :: VR
    Complex( wp ), Dimension(:,:), Intent(   Out ) :: VL

    Complex( wp ), Dimension( : ), Allocatable :: work

    Real( wp ), Dimension( : ), Allocatable :: rwork
    
    Integer :: d
    Integer :: INFO,LWORK

    d = Size( h, Dim = 1 )

    LWORK = 2*d
    Allocate( work( 1:lwork ) )
    Allocate( rwork( 1:2 * d ) )

    Call ZGEEV( 'V', 'V', d, h, d, W, VL, d, VR, d, WORK,LWORK, RWORK, INFO )
    Write( *, * ) 'Infor returned from diag ', info

  End Subroutine compl_diag_lr

End Module matrix_routines

Program NonHerm

  Use, Intrinsic :: iso_fortran_env, Only : wp => real64

  Use matrix_routines, Only : ham, compl_diag_lr

  Implicit None

  Integer,Parameter :: m = 20

  Integer :: i,j, k

  Complex( wp ), Allocatable :: h(:,:), W(:), evl(:,:), evr(:,:)

  Complex( wp ) :: t, g, h2(m,m), norm

  Allocate(h(m,m))
  Allocate(W(m))
  Allocate(evl(m,m))
  Allocate(evr(m,m))

  t = 1.0_wp;
  g = 0.0_wp

  Call ham(h,t,g)

  Call compl_diag_lr(h,W,evr,evl)

  Do j = 1, m
     norm = 0.0_wp
     Do k = 1, m
        norm = norm + Conjg(evl(k,j))*evr(k,j)
     End Do
     evr(:,j) = evr(:,j)/norm
  End Do

  h2 = Matmul(Conjg(Transpose(evl)), evr)
  Do i = 1, m
     h2( i, i ) = h2( i, i ) - 1.0_wp
  End Do

  Write( *, * ) 'The evals are'
  Do i = 1, m
     Write( *, '( i2, 4x,"( ", f16.12, ", ", f16.12, " )"  )' ) i, w( i )
  End Do

  Do j = 1, m
     Do i = 1, m
        If( Abs( h2( i, j ) ) > 1.0e-14 ) Then
           Write( *, '( "Non zero couple at ", i2, 1x, i2, &
                &". The Corresponding evals are ", 2( "( ", f12.8, ", ", f12.8, " )" : " and " ) )' ) &
                i, j, w( i ), w( j )
        End If
     End Do
  End Do

  Deallocate(h)
  Deallocate(W)
  Deallocate(evl)
  Deallocate(evr) 

End Program NonHerm

ijb@ijb-Latitude-5410:~/work/stack$ gfortran --version
GNU Fortran (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0
Copyright (C) 2019 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.

ijb@ijb-Latitude-5410:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -O -std=f2018 -Werror -Wuse-without-only -finit-real=snan -g zgeev.f90  -llapack
ijb@ijb-Latitude-5410:~/work/stack$ ./a.out
 Infor returned from diag            0
 The evals are
 1    (   2.000000000000,   0.000000000000 )
 2    (   1.902113032590,   0.000000000000 )
 3    (   1.618033988750,   0.000000000000 )
 4    (  -2.000000000000,   0.000000000000 )
 5    (  -1.902113032590,   0.000000000000 )
 6    (   1.902113032590,   0.000000000000 )
 7    (   1.618033988750,   0.000000000000 )
 8    (   1.175570504585,   0.000000000000 )
 9    (   1.175570504585,   0.000000000000 )
10    (   0.618033988750,   0.000000000000 )
11    (   0.618033988750,   0.000000000000 )
12    (  -0.000000000000,   0.000000000000 )
13    (   0.000000000000,   0.000000000000 )
14    (  -0.618033988750,   0.000000000000 )
15    (  -0.618033988750,   0.000000000000 )
16    (  -1.902113032590,   0.000000000000 )
17    (  -1.618033988750,   0.000000000000 )
18    (  -1.618033988750,   0.000000000000 )
19    (  -1.175570504585,   0.000000000000 )
20    (  -1.175570504585,   0.000000000000 )
Non zero couple at  3  7. The Corresponding evals are (   1.61803399,   0.00000000 ) and (   1.61803399,   0.00000000 )
Non zero couple at  8  9. The Corresponding evals are (   1.17557050,   0.00000000 ) and (   1.17557050,   0.00000000 )
ijb@ijb-Latitude-5410:~/work/stack$

Is biorthogonalization of left and right eigenvectors using ZGEEV in lapack guaranteed?

Answers (1)

Related Questions