python – zgeev()LAPACK的结果不正确/不一致
作者：互联网
我试图使用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和双精度,显式write()并将所有内容初始化为零以确保安全：
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]]
解决方法:
在程序中你有作为标量的rwork,它应该是一个大小为2 * N的数组,根据文档
更正此问题可以解决问题
标签：gfortran,python,lapack,fortran 来源： https://codeday.me/bug/20191003/1848052.html