Incorrect/inconsistent results from zgeev() LAPACK

Question

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

I assume you checked the values of `hamiltonian` (but it is not mentioned in the question). — albert, Jul 12 '18 at 10:42
Yes, I have included the print outs explicitly now. The Hamiltonian entering the functions each time are identical, though it does appear it is somehow being overwritten after calling zgeev() which I don't understand. — Chris, Jul 12 '18 at 10:48
The overwriting is explainable, from the documentation: `On exit, A has been overwritten.` — albert, Jul 12 '18 at 10:59
What are the `hamiltonian` values in python? Also check the values of the other arrays (and when not set do explicitly set the values to 0). — albert, Jul 12 '18 at 11:02
I've initialized everything as zero now to be safe, and printed info which is 0 as expected. Not sure what significance work has, though I've printed it also. I've also printed the Hamiltonian from Python, which is the same except for the column/row definition which I don't believe should make any difference. — Chris, Jul 12 '18 at 11:16

score 3 · Accepted Answer · answered Jul 12 '18 at 15:42

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

score 2 · Answer 2 · answered Jul 12 '18 at 10:21

2

According to the documentation of zgeev (on 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 
    )

The arrays you provide are of type complex(kind = 8) and not of type complex*16 and for RWord provided type real(kind = 8) instead of double precision. (Depending on the compiler the kind=8 can have a different meaning)

answered Jul 12 '18 at 10:21

albert

8,285
3
19
32

Thanks for your help. I was under the impression that for gfortran, though not guaranteed for other compilers, that complex(kind=8) and complex*16 were the same? Also I have just tested this and the issue remains (see edit). – Chris Jul 12 '18 at 10:28
1

They should be the same. Kind=8 might even be marginally better, because it is at least Fortran even if ugly and non-portable. The best thing is of course to avoid magic numbers and use named constants with some well-chosen value. – Vladimir F Героям слава Jul 12 '18 at 11:18
Though what you are doing is totally non-portable (You should never, NEVER, use magic numbers for kinds) they are almost certainly the same and the answer is not right. The kind of a complex number id the same as the kind of the real numbers that comprise it. Please see https://stackoverflow.com/questions/32935743/why-change-complex16-to-complex16-cause-the-runtime-increased-unreasonabl/32935983#32935983 – Ian Bush Jul 12 '18 at 11:19
@VladimirF In this case they might be identical but you are implicitly saying that the documentation / implementation of Lapack could be improved by using a better definition based on `selected_real_kind` (and friends). I agree with this. – albert Jul 12 '18 at 11:25
Are you suggesting this as a solution to the problem, or simply that kind=8 is bad form due to non-portability? I'm using Fortran with F2PY so thought kind=8 was safest, but I could switch over to using some form of integer, parameter :: r15 = selected_real_kind(15) and complex(kind=r15) if this is preferable – Chris Jul 12 '18 at 11:34
If your choice for kind parameter is constrained by **domain**, you should select it through `selected_real_kind`; if your choice is constrained by **storage**, you should select using the constants and procedures in `iso_fortran_env`, like `real64` or `nuneric_storage_size`. – Rodrigo Rodrigues Jul 12 '18 at 12:59
Thanks, I'll bear that in mind for the future. However, it doesn't solve my current problem of zgeev() giving incorrect results. – Chris Jul 12 '18 at 13:42
Isn't it necessary for Rwork to be an array..? (it seems a scaler in OP's code) – roygvib Jul 12 '18 at 15:35
I think this might solve the problem, from the documentation: `RWORK is DOUBLE PRECISION array, dimension (2*N)` – albert Jul 12 '18 at 15:41

Incorrect/inconsistent results from zgeev() LAPACK

2 Answers2

Linked

Related