Resultados incorrectos / inconsistentes de zgeev () LAPACK

Estoy intentando usar ZGEEV para calcular los valores propios y los vectores propios, sin embargo, estoy teniendo algunos problemas con la salida como incorrecta y también inconsistente cuando se usa a diferentes niveles de optimización. A continuación se muestra mi código Fortran con resultados en los niveles de optimización de -O1 y -O2. También he incluido el código de Python para la comparación.

Solo puedo asumir que estoy llamando a zgeev() incorrectamente de alguna manera, sin embargo no puedo determinar cómo. Creo que es poco probable que sea un problema con mi instalación de LAPACK, ya que he comparado la salida en dos computadoras diferentes, en Windows y Linux.

Código de 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 

Resultados en -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) 

Resultados en -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) 

    Código 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) 

    Resultados de Python:

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

    EDITAR:

    Usando complex * 16 y doble precisión como se especifica en la documentación, escriba () explícitamente e inicie todo como cero para estar seguro:

     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 

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

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

    Pitón:

     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) 

    Resultado:

     hamiltonian [[ 1000000. 250000.] [ 750000. -1000000.]] hamiltonian [[ 1000000. 250000.] [ 750000. -1000000.]] eigvals [ 1089724.73588517 -1089724.73588517] eigvectors [[ 0.94121724 -0.11878597] [ 0.33780187 0.99291988]] 

    En el progtwig que tiene que trabajar como un escalar, debe ser una matriz de tamaño 2 * N de acuerdo con la documentación en

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

    Corrigiendo esto arregla el problema.

    De acuerdo con la documentación de zgeev (en http://www.netlib.org/lapack/explore-html/db/d55/group__complex16_g_eeigen_ga0eb4e3d75621a1ce1685064d64c1ac58f0.html#ga0eb4e3d75621a1ce1686464641

     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 ) 

    Los arreglos que proporcione son de tipo complex(kind = 8) y no de tipo complex*16 y para RWord proporcionado tipo real(kind = 8) lugar de double precision . (Dependiendo del comstackdor, el kind=8 puede tener un significado diferente)