Lapack, DGEEV does not give me the right eigenvalues

When trying to calculate the eigenvalues of a matrix with the DGEEV function, I noticed that they were non conjugate and thus wrong.

I tried testing with the simple (0 1; -1 0) matrix (which has eigenvalues i and -i), but it doesn't give me the right solution.

Can someone point me to the error in my code?

program myProgram
    real(kind=dp), dimension(2) :: ew_real, ew_imag, ev
    real(kind=dp), dimension(4*5,4*5) :: work 
    integer info
    real(kind=dp), DIMENSION(2, 2) :: array
    array = reshape((/ 0.0, -1.0, 1.0, 0.0/), shape(array))      

    call DGEEV('N','V',2,array,2,ew_real,ew_imag,ev,2,ev,2,work,size(work,1),info)
end program

Output of my program: (0.00000,0.00000) and (0.00000,-0.40824)

Solution

Mistakes

The ev array is only of dimension(2) whereas it should be of dimension(2,2).
work should be a one-dimensional array of length 4n=4*2.

Note that the first one actually corrupts your memory and thus, you get the wrong eigenvalues.

Improvement

Moreover, you should be using the lapack95 module which gives generic interfaces that are easier to call.

The following example shows both calls

program main
  use iso_fortran_env, only: dp => real64
  use lapack95

  integer                       :: info
  real(kind=dp), dimension(2)   :: ew_real, ew_imag
  real(kind=dp), dimension(4*2) :: work
  real(kind=dp), dimension(2,2) :: array, ev_left, ev_right
  array = reshape([0, -1, 1, 0], shape(array))

  ! lapack
  call dgeev('N','V',2,array,2,ew_real,ew_imag,ev_left,2,ev_right,2,work,size(work),info)
  print *, 'eval1', cmplx(ew_real(1), y=ew_imag(1))
  print *, 'eval2', cmplx(ew_real(2), y=ew_imag(2))

  ! lapack95 call
  call geev(array, ew_real, ew_imag, vr=ev_right)
  print *, 'eval1', cmplx(ew_real(1), y=ew_imag(1))
  print *, 'eval2', cmplx(ew_real(2), y=ew_imag(2))
end program

The output is

$ ./main
 eval1               (0.0000000000000000,1.0000000000000000)
 eval2              (0.0000000000000000,-1.0000000000000000)
 eval1               (0.0000000000000000,1.0000000000000000)
 eval2              (0.0000000000000000,-1.0000000000000000)