重复时 LAPACK 中 ZGEEV 的特征值不正确
Incorrect eigenvalues with ZGEEV in LAPACK when repeated
我有以下代码可以找到非 Hermitian 复数矩阵的特征值 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
在这段代码中,我使用了一个耗散或不对称参数 lambda
,我将其固定为 l=2.d0
和 运行 一个 do-loop 三次。我应该在每次迭代中获得相同的特征值。令我惊讶的是,前两个循环给出了相同的特征值集,但它们的顺序不同,第三个循环生成了一组完全不同的特征值!
可能的错误是什么?
输出:
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.: 有一个 有类似的反对意见,但我不认为我犯了错误定义 RWORK
的错误 RWORK
=28=]。所以请不要立即将我的问题标记为重复。
问题是,引用 Zgeev 页面
A is COMPLEX*16 array, dimension (LDA,N)
On entry, the N-by-N matrix A.
On exit, A has been overwritten.
重点是我的。因此,您需要每次都重新初始化整个矩阵,而不仅仅是非零位。修复此问题,使您的代码符合标准,并将其移动到更符合现代实践的样式我得到
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$
我有以下代码可以找到非 Hermitian 复数矩阵的特征值 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
在这段代码中,我使用了一个耗散或不对称参数 lambda
,我将其固定为 l=2.d0
和 运行 一个 do-loop 三次。我应该在每次迭代中获得相同的特征值。令我惊讶的是,前两个循环给出了相同的特征值集,但它们的顺序不同,第三个循环生成了一组完全不同的特征值!
可能的错误是什么?
输出:
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.: 有一个 RWORK
的错误 RWORK
=28=]。所以请不要立即将我的问题标记为重复。
问题是,引用 Zgeev 页面
A is COMPLEX*16 array, dimension (LDA,N) On entry, the N-by-N matrix A. On exit, A has been overwritten.
重点是我的。因此,您需要每次都重新初始化整个矩阵,而不仅仅是非零位。修复此问题,使您的代码符合标准,并将其移动到更符合现代实践的样式我得到
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$