title: Incorrect/inconsistent results from zgeev() LAPACK
body: <p>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. </p>
  
  <p>I can only assume that I am calling <code>zgeev()</code> 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.</p>
  
  <p>Fortran code:</p>
  
  <pre><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
  </code></pre>
  
  <p>Results at -O1:</p>
  
  <pre><code> 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)
  </code></pre>
  
  <p>Results at -O2:</p>
  
  <pre><code> 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)
  </code></pre>
  
  <p>Python code:</p>
  
  <pre><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)
  </code></pre>
  
  <p>Python results:</p>
  
  <pre><code>eigvals [ 1089724.73588517 -1089724.73588517]
  eigvectors [[ 0.94121724 -0.11878597] [ 0.33780187  0.99291988]]
  </code></pre>
  
  <p><strong>EDIT:</strong></p>
  
  <p>Using complex*16 and double precision as specified in documentation, explicit write() and initializing everything as zero to be safe:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>Output -O1:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>Output -O2:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>Python:</p>
  
  <pre><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)
  
  print('hamiltonian', hamiltonian)
  
  eigvals, eigvectors = np.linalg.eig(hamiltonian)
  
  print('hamiltonian', hamiltonian)
  print('eigvals', eigvals)
  print('eigvectors', eigvectors)
  </code></pre>
  
  <p>Result:</p>
  
  <pre><code>hamiltonian [[ 1000000.   250000.] [  750000. -1000000.]]
  hamiltonian [[ 1000000.   250000.] [  750000. -1000000.]]
  eigvals [ 1089724.73588517 -1089724.73588517]
  eigvectors [[ 0.94121724 -0.11878597] [ 0.33780187  0.99291988]]
  </code></pre>