Incorrect/inconsistent 结果来自 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 和双精度、显式 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]]

Answer 1

根据 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 可能有不同的含义）

Answer 2

在您将 rwork 作为标量的程序中，根据

上的文档，它应该是一个大小为 2*N 的数组

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

更正此问题即可解决

Incorrect/inconsistent 结果来自 zgeev() LAPACK

Incorrect/inconsistent results from zgeev() LAPACK

python

fortran

gfortran

lapack