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

More details please - in particular a *complete*, *compilable* code that illustrates the problem you are having, sample output, explanation why that output is not what you expect, and an explanation of hat you do expect. — Ian Bush, May 26 '22 at 18:18
I have added a code. Please look into it. You can tune the Non-Hermitianity through the parameter g. — Nil, May 26 '22 at 21:45

Ian Bush · Accepted Answer · 2022-05-27T07:30:04.753

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$

If the parameter 'g' is non-zero, then the degeneracy is broken. Are you saying that biorthogonality is still not guaranteed by LAPACK between the two sets of eigenvectors? If you know of any alternative that would help. — Nil, May 27 '22 at 07:27
Changing g does *not* break the degeneracy - it simply shifts the eigenvalue spectrum uniformly. Changing the diagonal to random numbers does break the degenracy, and in that case biorthogonality is strictly observed. Try it! — Ian Bush, May 27 '22 at 07:33
If you really need to enforce biorthog even in the degenerate case 2 seconds of searching shows up https://joshuagoings.com/2015/04/03/biorthogonalizing-left-and-right-eigenvectors-the-easy-lazy-way/ — Ian Bush, May 27 '22 at 07:36
Note for large matrices this will not be especially efficient, but it will work. Really you only need to orthogonalise in the space containing the "problem" evecs, rather than the whole space as this does -and the "problem" space will usually be much smaller. That said you've just done a diag, so a matmul and LU decomp will not be a huge overhead. — Ian Bush, May 27 '22 at 08:17

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

1 Answers1