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

2022-01-14 00:00:00 python fortran gfortran lapack

问题描述

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

I am attempting to use ZGEEV to calculate eigenvalues and eigenvectors, however am having some trouble with the output being incorrect and also inconsistent when used at different optimization levels. Below is my Fortran code with results at -O1 and -O2 optimization levels. I have also included Python code for comparison.

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

I can only assume that I am calling zgeev() incorrectly somehow, however I am not able to determine how. I believe it is unlikely to be an issue with my LAPACK installation as I have compared the output on two different computers, on Windows and Linux.

Fortran 代码:

Fortran code:

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() 并将所有内容初始化为零以确保安全:

Using complex*16 and double precision as specified in documentation, explicit write() and initializing everything as zero to be safe:

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 的数组

In the program you have rwork as a scalar, it should be an array of size 2*N according to the documentation at

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

纠正这个问题可以解决问题

Correcting this fixes the problem

相关文章