dgeev_rel Subroutine

   pure subroutine dgeev_rel(matrix, eig_vecr, eig_val, eig_vecl)
      real(rk), dimension(:,:), intent(in)                         :: matrix
      real(rk), dimension(:),   allocatable, intent(out)           :: eig_val
      real(rk), dimension(:,:), allocatable, intent(out)           :: eig_vecr
      real(rk), dimension(:,:), allocatable, intent(out), optional :: eig_vecl
      real(rk), dimension(size(matrix,1))                          :: wr, wi
      real(rk), dimension(size(matrix,1),size(matrix,1))           :: vl, vr
      real(rk), dimension(size(matrix,1),size(matrix,1))           :: A
      real(rk), dimension(:), allocatable                          :: work
      integer                                                      :: m, lwork, info
      real(rk)                                                     :: work1(1)

      interface
         pure subroutine dgeev(fjobvl, fjobvr, fn, fA, flda, fwr, fwi, fvl, fldvl, fvr, fldvr, fwork, flwork, finfo)
            import :: rk
            character, intent(in)    :: fjobvl, fjobvr
            integer,   intent(in)    :: fn, flda, fldvl, fldvr, flwork, finfo
            real(rk),  intent(inout) :: fA(flda, *)
            real(rk),  intent(out)   :: fwr(fn), fwi(fn), fvl(fldvl, *), fvr(fldvr, *), fwork(flwork)
         end subroutine
      end interface

      m = size(matrix, 1)
      A = matrix

      if (present(eig_vecl)) then
         eig_vecl = matrix
         call dgeev('V', 'V', m, A, m, wr, wi, vl, m, vr, m, work1, -1, info)
      else
         call dgeev('N', 'V', m, A, m, wr, wi, vl, m, vr, m, work1, -1, info)
      end if

      lwork = nint(work1(1))
      allocate(work(lwork))

      if (present(eig_vecl)) then
         call dgeev('V', 'V', m, A, m, wr, wi, vl, m, vr, m, work, lwork, info)
         eig_vecl = vl
      else
         call dgeev('N', 'V', m, A, m, wr, wi, vl, m, vr, m, work, lwork, info)
      end if

      eig_val  = wr
      eig_vecr = vr

      deallocate(work)

   end subroutine dgeev_rel