zgeev()LAPACK的结果不正确/不一致

Question

我正试图使用ZGEEV来计算特征值和特征向量，但在不同的优化级别使用时，输出不正确，也不一致，因此遇到了一些问题。下面是我的Fortran代码，其中包含-O1和-O2优化级别的结果。我还包含了用于比较的Python代码。

我只能假设我以某种方式错误地调用了zgeev()，但我无法确定如何调用。我相信这不太可能是我的LAPACK安装的问题，因为我已经比较了Windows和Linux上两台不同计算机的输出。

Fortran代码：

program example_main
use example_subroutine
implicit none
complex(kind = 8) :: eigval(2), dummy(2, 2), work(4), eig_vector(2, 2)
real(kind = 8) :: Rwork
complex(kind = 8), dimension(2, 2) :: hamiltonian
integer :: info, count
call calculate_hamiltonian(hamiltonian)
call ZGEEV('N', 'V', 2, hamiltonian, 2, eigval, dummy, 4, eig_vector, 2, work, 4, Rwork, info)
end program example_main
module example_subroutine
contains
subroutine calculate_hamiltonian(hamiltonian)
implicit none
integer :: count
complex(kind = 8), dimension(2, 2), intent(out) :: hamiltonian
complex(kind = 8), dimension(2, 2) :: spin_x, spin_z
spin_x = 0.5 * (reshape((/ 0.D0, 1.D0, 1.D0, 0.D0/), shape(spin_x), order = (/2, 1/)))
spin_z = 0.5 * (reshape((/ 1.D0, 0.D0, 0.D0, -1.D0/), shape(spin_z), order = (/2, 1/)))
hamiltonian =  2D6 * spin_z + 1D6 * spin_x + 1E6 * matmul(spin_x, spin_z)
end subroutine calculate_hamiltonian
end module

-O1:的结果

eigval               
(1089724.7358851689,0.0000000000000000) (-1089724.7358851684,0.0000000000000000)
eig_vector               
(1.0000000000000000,0.0000000000000000) (0.0000000000000000,-0.0000000000000000) (1.0000000000000000,0.0000000000000000) (0.0000000000000000,0.0000000000000000)

O2下的结果：

eigval         
(1089724.7358851689,1.20522527882675885E-014) (0.99999999999998823,0.0000000000000000)
eig_vector        
(2.55688391396797063E-006,-0.0000000000000000) (0.99999999999673128,0.0000000000000000) (-1.09782752690336509E-007,0.0000000000000000) (0.99999999999999412,0.0000000000000000)

Python代码：

spin_x = 1/2 * np.array([[0, 1],  [1, 0]])
spin_z = 1/2 * np.array([[1, 0], [0, -1]])
hamiltonian = 2E6 * spin_z + 1E6 * spin_x + 1E6 * np.matmul(spin_x, spin_z)
eigvals, eigvectors = np.linalg.eig(hamiltonian)

Python结果：

eigvals [ 1089724.73588517 -1089724.73588517]
eigvectors [[ 0.94121724 -0.11878597] [ 0.33780187  0.99291988]]

编辑：

使用文档中指定的复数*16和双精度，显式写入((并将所有内容初始化为零以确保安全：

module example_subroutine
contains
subroutine calculate_hamiltonian(hamiltonian)
implicit none
complex*16, dimension(2, 2), intent(out) :: hamiltonian
complex*16, dimension(2, 2) :: spin_x, spin_z
hamiltonian = 0
spin_x = 0
spin_z = 0
spin_x = 0.5 * (reshape((/ 0.D0, 1.D0, 1.D0, 0.D0/), shape(spin_x), order = (/2, 1/)))
spin_z = 0.5 * (reshape((/ 1.D0, 0.D0, 0.D0, -1.D0/), shape(spin_z), order = (/2, 1/)))
hamiltonian =  2D6 * spin_z + 1D6 * spin_x + 1E6 * matmul(spin_x, spin_z)
write(6, *) 'hamiltonian', hamiltonian
end subroutine calculate_hamiltonian
end module
program example_main
use example_subroutine
implicit none
complex*16 :: eigval(2), dummy(2, 2), work(4), eig_vector(2, 2)
double precision :: Rwork
complex*16, dimension(2, 2) :: hamiltonian
integer :: info
eigval = 0
dummy = 0
work = 0
eig_vector = 0
Rwork = 0
info = 0
hamiltonian = 0
call calculate_hamiltonian(hamiltonian)
write(6, *) 'hamiltonian before', hamiltonian
call ZGEEV('N', 'V', 2, hamiltonian, 2, eigval, dummy, 4, eig_vector, 2, work, 4, Rwork, info)
write(6, *) 'hamiltonian after', hamiltonian
write(6, *) 'eigval', eigval
write(6, *) 'eig_vector', eig_vector
write(6, *) 'info', info
write(6, *) 'work', work
end program example_main

输出-O1:

hamiltonian    
(1000000.0000000000,0.0000000000000000) (750000.00000000000,0.0000000000000000) (250000.00000000000,0.0000000000000000) (-1000000.0000000000,0.0000000000000000)
hamiltonian before      
(1000000.0000000000,0.0000000000000000) (750000.00000000000,0.0000000000000000) (250000.00000000000,0.0000000000000000) (-1000000.0000000000,0.0000000000000000)
hamiltonian after          
(0.99999999999999989,0.0000000000000000) (0.0000000000000000,0.0000000000000000) (500000.00000000012,0.0000000000000000) (-1089724.7358851684,0.0000000000000000)
eigval               
(1089724.7358851689,0.0000000000000000) (-1089724.7358851684,0.0000000000000000)
eig_vector        
(1.0000000000000000,0.0000000000000000) (0.0000000000000000,-0.0000000000000000) (1.0000000000000000,0.0000000000000000) (0.0000000000000000,0.0000000000000000)
info           0
work               
(260.00000000000000,0.0000000000000000) (-1089724.7358851684,0.0000000000000000) (1.0000000000000000,0.0000000000000000) (1.0000000000000000,0.0000000000000000)

输出-O2:

hamiltonian               
(1000000.0000000000,0.0000000000000000) (750000.00000000000,0.0000000000000000) (250000.00000000000,0.0000000000000000) (-1000000.0000000000,0.0000000000000000)
hamiltonian before               
(1000000.0000000000,0.0000000000000000) (750000.00000000000,0.0000000000000000) (250000.00000000000,0.0000000000000000) (-1000000.0000000000,0.0000000000000000)
hamiltonian after               
(1089724.7358851689,0.0000000000000000) (0.0000000000000000,0.0000000000000000) (500000.00000000012,0.0000000000000000) (-1089724.7358851684,0.0000000000000000)
eigval         
(1089724.7358851689,1.20522527882675885E-014) (0.99999999999998823,0.0000000000000000)
eig_vector        
(2.55688391396797063E-006,-0.0000000000000000) (0.99999999999673128,0.0000000000000000) (-1.09782752690336509E-007,0.0000000000000000) (0.99999999999999412,0.0000000000000000)
info           0
work               
(260.00000000000000,0.0000000000000000) (-1089724.7358851684,0.0000000000000000) (1.0000000000000000,0.0000000000000000) (1.0000000000000000,0.0000000000000000)

Python:

spin_x = 1/2 * np.array([[0, 1],  [1, 0]])
spin_z = 1/2 * np.array([[1, 0], [0, -1]])
hamiltonian = 2E6 * spin_z + 1E6 * spin_x + 1E6 * np.matmul(spin_x, spin_z)
print('hamiltonian', hamiltonian)
eigvals, eigvectors = np.linalg.eig(hamiltonian)
print('hamiltonian', hamiltonian)
print('eigvals', eigvals)
print('eigvectors', eigvectors)

结果：

hamiltonian [[ 1000000.   250000.] [  750000. -1000000.]]
hamiltonian [[ 1000000.   250000.] [  750000. -1000000.]]
eigvals [ 1089724.73588517 -1089724.73588517]
eigvectors [[ 0.94121724 -0.11878597] [ 0.33780187  0.99291988]]

Answer 1

根据的文档，在您将rwork作为标量的程序中，它应该是一个大小为2*N的数组

http://www.netlib.org/lapack/explore-html/db/d55/group__complex16_g_eeigen_ga0eb4e3d75621a1ce1685064db1ac58f0.html#ga0eb4e3d75621a1ce1685064db1ac58f0

更正此问题可修复问题

Answer 2

根据zgeev(http://www.netlib.org/lapack/explore-html/db/d55/group__complex16_g_eeigen_ga0eb4e3d75621a1ce1685064db1ac58f0.html#ga0eb4e3d75621a1ce1685064db1ac58f0)：

subroutine zgeev    (   character   JOBVL,
character   JOBVR,
integer     N,
complex*16, dimension( lda, * )     A,
integer     LDA,
complex*16, dimension( * )      W,
complex*16, dimension( ldvl, * )    VL,
integer     LDVL,
complex*16, dimension( ldvr, * )    VR,
integer     LDVR,
complex*16, dimension( * )      WORK,
integer     LWORK,
double precision, dimension( * )    RWORK,
integer     INFO 
)

您提供的数组的类型为complex(kind = 8)，而不是complex*16，并且为RWord提供的数组类型为real(kind = 8)，而不是double precision。(根据编译器的不同，kind=8可能有不同的含义(