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
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$