是否保证在 lapack 中使用 ZGEEV 对左右特征向量进行双正交化？
Is biorthogonalization of left and right eigenvectors using ZGEEV in lapack guaranteed?
假设我想对角化 H = H0 + iV 形式的复杂非厄米矩阵，
其中 H0 和 V 是 Hermitian 矩阵。设R和L是分别包含右特征向量和左特征向量作为列的矩阵。则根据双正交化条件应满足 conjugate(transpose(L)).R = I。但是使用 ZGEEV 得到的 L 和 R 似乎不满足这个条件。
请查看以下 fortran90 代码：
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
您的矩阵具有退化特征值。对应于此类 eval 的 evec 之间不保证双正交性。就是这些情况non-zero，所以LAPACK没有错。注意在这种情况下可以强制执行双正交性，但是 LAPACK documentation 确实 而不是 说它会做这样的事情，所以你不能依赖它 - 正如你在这里看到的那样。如果您需要它，您将不得不编写自己的代码来完成它。
另请注意，您的代码不是标准 Fortran - complex*16 和类似代码不是也从来不是 Fortran 的一部分，不应使用。请了解 Fortran kinds 。还应弃用外部子例程，改用包含的或最好的模块子程序。这是您的程序的标准符合版本，具有更现代的编程风格，并且更易于阅读输出：
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$
假设我想对角化 H = H0 + iV 形式的复杂非厄米矩阵， 其中 H0 和 V 是 Hermitian 矩阵。设R和L是分别包含右特征向量和左特征向量作为列的矩阵。则根据双正交化条件应满足 conjugate(transpose(L)).R = I。但是使用 ZGEEV 得到的 L 和 R 似乎不满足这个条件。 请查看以下 fortran90 代码：
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
您的矩阵具有退化特征值。对应于此类 eval 的 evec 之间不保证双正交性。就是这些情况non-zero，所以LAPACK没有错。注意在这种情况下可以强制执行双正交性，但是 LAPACK documentation 确实 而不是 说它会做这样的事情，所以你不能依赖它 - 正如你在这里看到的那样。如果您需要它，您将不得不编写自己的代码来完成它。
另请注意，您的代码不是标准 Fortran - complex*16 和类似代码不是也从来不是 Fortran 的一部分，不应使用。请了解 Fortran kinds 。还应弃用外部子例程，改用包含的或最好的模块子程序。这是您的程序的标准符合版本，具有更现代的编程风格，并且更易于阅读输出：
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$