! test_eigen2.f90  use the LAPACK eigenvalue code  zgeev.f to
!                 test various eigenvalue definitions/identies
!                 Av = ev   det|A - eI| =0  
!                 B = P A P^-1  A and B have same eigenvalues
!                 ramification A with row and/or column interchange, same e
!                 cA has eigenvalues ce
!                 largest e <= L1 norm (max of sum of abs of values in row or col)
!                 real positive definate matrix A has all real e
!                 ramification real symmetric A has all real e

program test_eigen2
  implicit none
  integer, parameter :: n = 4
  integer, parameter :: m = 4
  double complex, dimension(n,m) :: A ! my matricies for testing
  double complex, dimension(n,m) :: B
  double complex, dimension(n,m) :: C
  double complex, dimension(n,m) :: D
  double complex, dimension(n)   :: e ! my vectors for testing
  double complex, dimension(n)   :: u
  double complex, dimension(n)   :: v
  double precision               :: L1norm
  double precision               :: paste
  integer                        :: i, j, k 

  ! stuff needed by  zgeev
  double complex, dimension(n,n)   :: vR
  double complex, dimension(n,n)   :: vL
  double complex, dimension(n*n)   :: WORK 
  double complex, dimension(n*n)   :: LWORK
  double precision, dimension(2*n) :: RWORK
  integer                          :: info
  
  interface
    function znorm1(n, m, A) result(L1norm)
      integer, intent(in) :: n, m
      double complex, dimension(n,m), intent(in) :: A
      double precision :: L1norm
    end function znorm1
    subroutine zvecprint(v, c)
      double complex, dimension(:), intent(in) :: v
      character(len=*), intent(in) :: c
    end subroutine zvecprint
  end interface 
  
  print *, "test_eigen2.f90"
  do i = 1,n
    do j = 1,m
      A(i,j) = dcmplx(max(i,j),max(5-i,5-j))
    end do
  end do
  call zmatprint(n, m, A, "A")
  L1norm = znorm1(n, m, A)
  print *, "L1norm=", L1norm
  call zgeev('V', 'V', 4, A,  4, e, vL,  4,  vR,  4, WORK, 16,  RWORK, info)
            ! VL   VR  N  A  LDA W  VL LDVL  VR LDVR WORK LWORK RWORK  INFO
  call zvecprint(e, "eigenvalues")
  call zmatprint(n, m, vR, "right values")
  call zmatprint(n, m, vL, "left values")
  call evalvec(n, n, A, e, vR)
  call evalvec(n, n, A, e, vL)
  call zcheckeigen(n, A, e, vR);
  call zcheckeigen(n, A, e, vL);

end program test_eigen2

! support routines 

subroutine zCheckEigen(n, A, e, V)
  implicit none
  integer, intent(in) :: n
  double complex, dimension(n,n), intent(in) :: A ! matrix
  double complex, dimension(n),   intent(in) :: e ! eigenvalues
  double complex, dimension(n,n), intent(in) :: V ! eigenvectors

  double precision :: error, max_error, mean_error
  double complex, dimension(n,n) :: B
  double complex, dimension(n,n) :: Ident
  double complex, dimension(n)   :: vk
  double complex, dimension(n)   :: ve
  double complex, dimension(n)   :: Av
  double complex, dimension(n)   :: Avve
  integer                        :: i, j, k, m

  interface
    subroutine zvecprint(v, c)
      double complex, dimension(:), intent(in) :: v
      character(len=*), intent(in) :: c
    end subroutine zvecprint
  end interface 

  print *, "zCheckEigen"
  do k = 1,n
    do j = 1,n
      vk(j) = V(j,k)
      ve(j) = e(k)*V(j,k)
    end do
    call zMulMatVec(n, n, A, vk, Av);
    call zvecprint(vk, "v ")
    call zvecprint(ve, "ev")
    call zvecprint(Av, "Av")
    do m =1,n
      Avve(m) = Av(m)-Ve(m)
    end do
    call zvecprint(Avve, "Err")
  end do
end subroutine zCheckEigen

subroutine zMulMatVec(n, m, A, v, u)
  implicit none
  integer, intent(in) :: n, m
  double complex, dimension(n,m), intent(in) :: A
  double complex, dimension(n), intent(in) :: v
  double complex, dimension(m), intent(inout) :: u
  integer :: i, j

  do i = 1,n
    u(i) = 0.0
    do j = 1,m
      u(i) = u(i) + A(i,j)*v(j)
    end do
  end do
end subroutine zMulMatVec

function znorm1(n, m, A) result(L1norm)
  implicit none
  integer, intent(in) :: n, m
  double complex, dimension(n,m), intent(in) :: A
  double precision :: L1norm
  double precision :: temp
  integer :: i, j
  L1norm = 0.0
  do i=1,n
    temp = 0.0
    do j=1,m
      temp = temp + abs(A(i,j))
    end do
    if (temp > L1norm) L1norm = temp
  end do
end function znorm1

subroutine zvecprint(v, c)
  implicit none
  double complex, dimension(:), intent(in) :: v
  character(len=*), intent(in) :: c
  integer :: i, n

  n = size(v)
  do i=1,n
    print *, c, "(", i, ") = ", v(i)
  end do
end subroutine zvecprint

subroutine zmatprint(n, m, A, c)
  implicit none
  integer, intent(in) :: n, m
  double complex, dimension(n,m), intent(in) :: A
  character, intent(in) :: c
  integer :: i, j
  do i=1,n
    do j=1,m
      print *, c, "(", i, ",", j, ") = ", A(i,j)
    end do
  end do
end subroutine zmatprint

subroutine cmatprint(n, m, A, c)
  implicit none
  integer, intent(in) :: n, m
  complex, dimension(n,m), intent(in) :: A
  character, intent(in) :: c
  integer :: i, j
  do i=1,n
    do j=1,m
      print *, c, "(", i, ",", j, ") = ", A(i,j)
    end do
  end do
end subroutine cmatprint

subroutine smatprint(n, m, A, c)
  implicit none
  integer, intent(in) :: n, m
  real, dimension(n,m), intent(in) :: A
  character, intent(in) :: c
  integer :: i, j
  do i=1,n
    do j=1,m
      print *, c, "(", i, ",", j, ") = ", A(i,j)
    end do
  end do
end subroutine smatprint

subroutine dmatprint(n, m, A, c)
  implicit none
  integer, intent(in) :: n, m
  double precision, dimension(n,m), intent(in) :: A
  character, intent(in) :: c
  integer :: i, j
  do i=1,n
    do j=1,m
      print *, c, "(", i, ",", j, ") = ", A(i,j)
    end do
  end do
end subroutine dmatprint
  
!     SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
!    $                  WORK, LWORK, RWORK, INFO )
!
!  -- LAPACK driver routine (version 3.0) --
!     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
!     Courant Institute, Argonne National Lab, and Rice University
!     June 30, 1999
!
!     .. Scalar Arguments ..
!     CHARACTER          JOBVL, JOBVR
!     INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
!     ..
!     .. Array Arguments ..
!     DOUBLE PRECISION   RWORK( * )
!     COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
!    $                   W( * ), WORK( * )
!     ..
!
!  Purpose
!  =======
!
!  ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
!  eigenvalues and, optionally, the left and/or right eigenvectors.
!
!  The right eigenvector v(j) of A satisfies
!                   A * v(j) = lambda(j) * v(j)
!  where lambda(j) is its eigenvalue.
!  The left eigenvector u(j) of A satisfies
!                u(j)**H * A = lambda(j) * u(j)**H
!  where u(j)**H denotes the conjugate transpose of u(j).
!
!  The computed eigenvectors are normalized to have Euclidean norm
!  equal to 1 and largest component real.
!
!  Arguments
!  =========
!
!  JOBVL   (input) CHARACTER*1
!          = 'N': left eigenvectors of A are not computed;
!          = 'V': left eigenvectors of are computed.
!
!  JOBVR   (input) CHARACTER*1
!          = 'N': right eigenvectors of A are not computed;
!          = 'V': right eigenvectors of A are computed.
!
!  N       (input) INTEGER
!          The order of the matrix A. N >= 0.
!
!  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
!          On entry, the N-by-N matrix A.
!          On exit, A has been overwritten.
!
!  LDA     (input) INTEGER
!          The leading dimension of the array A.  LDA >= max(1,N).
!
!  W       (output) COMPLEX*16 array, dimension (N)
!          W contains the computed eigenvalues.
!
!  VL      (output) COMPLEX*16 array, dimension (LDVL,N)
!          If JOBVL = 'V', the left eigenvectors u(j) are stored one
!          after another in the columns of VL, in the same order
!          as their eigenvalues.
!          If JOBVL = 'N', VL is not referenced.
!          u(j) = VL(:,j), the j-th column of VL.
!
!  LDVL    (input) INTEGER
!          The leading dimension of the array VL.  LDVL >= 1; if
!          JOBVL = 'V', LDVL >= N.
!
!  VR      (output) COMPLEX*16 array, dimension (LDVR,N)
!          If JOBVR = 'V', the right eigenvectors v(j) are stored one
!          after another in the columns of VR, in the same order
!          as their eigenvalues.
!          If JOBVR = 'N', VR is not referenced.
!          v(j) = VR(:,j), the j-th column of VR.
!
!  LDVR    (input) INTEGER
!          The leading dimension of the array VR.  LDVR >= 1; if
!          JOBVR = 'V', LDVR >= N.
!
!  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
!          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!
!  LWORK   (input) INTEGER
!          The dimension of the array WORK.  LWORK >= max(1,2*N).
!          For good performance, LWORK must generally be larger.
!
!          If LWORK = -1, then a workspace query is assumed; the routine
!          only calculates the optimal size of the WORK array, returns
!          this value as the first entry of the WORK array, and no error
!          message related to LWORK is issued by XERBLA.
!
!  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
!
!  INFO    (output) INTEGER
!          = 0:  successful exit
!          < 0:  if INFO = -i, the i-th argument had an illegal value.
!          > 0:  if INFO = i, the QR algorithm failed to compute all the
!                eigenvalues, and no eigenvectors have been computed;
!                elements and i+1:N of W contain eigenvalues which have
!                converged.
!