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]]
In the program you have rwork as a scalar, it should be an array of size 2*N according to the documentation at
Correcting this fixes the problem