Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
246 views
in Technique[技术] by (71.8m points)

python - Incorrect/inconsistent results from zgeev() LAPACK

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.

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 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

Results at -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)

Results at -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 code:

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 results:

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

EDIT:

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

Output -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)

Output -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)

Result:

hamiltonian [[ 1000000.   250000.] [  750000. -1000000.]]
hamiltonian [[ 1000000.   250000.] [  750000. -1000000.]]
eigvals [ 1089724.73588517 -1089724.73588517]
eigvectors [[ 0.94121724 -0.11878597] [ 0.33780187  0.99291988]]
See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

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


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...