我正在尝试使用ZGEEV计算特征值和特征向量,但是在以不同的优化级别使用时,输出不正确且不一致也有些麻烦。以下是我的Fortran代码,结果为-O1和-O2优化级别。我还提供了Python代码进行比较。
我只能假设我打错了电话zgeev(),但是我无法确定如何打电话。我相信我的LAPACK安装不太可能出现问题,因为我已经比较了Windows和Linux上两台不同计算机上的输出。
zgeev()
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]]
编辑:
使用文档中指定的complex * 16和double precision,显式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)
蟒蛇:
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的数组
http://www.netlib.org/lapack/explore- html/db/d55/group__complex16_g_eeigen_ga0eb4e3d75621a1ce1685064db1ac58f0.html#ga0eb4e3d75621a1ce1685064db1ac58f0
更正此问题即可解决