! test_eigen.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_eigen
  implicit none
  integer :: n=4, m=4
  double complex, dimension(4,4) :: A ! my matricies for testing
  double complex, dimension(4,4) :: B
  double complex, dimension(4,4) :: C
  double complex, dimension(4,4) :: D
  double complex, dimension(4)   :: e ! my vectors for testing
  double complex, dimension(4)   :: u
  double complex, dimension(4)   :: v
  double precision               :: L1norm
  double precision               :: paste

  ! stuff needed by  zgeev
  double complex, dimension(4,4) :: vR
  double complex, dimension(4,4) :: vL
  double complex, dimension(16)  :: WORK
  double complex, dimension(16)  :: LWORK
  double precision, dimension(8) :: 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
  end interface

  print *, "test_eigen.f90"
  A = reshape((/ 85, 71, 66, 44,  76, 45, 45, 62,    &   
                 90, 50, 21, 17,  60, 80, 55, 28 /), shape(A))
  call zmatprint(n, m, A, "A")
  L1norm = znorm1(n, m, A)
  print *, "L1norm=", L1norm
  v = reshape((/(1.0,2.0), (3.0,4.0), (5.0,6.0), (7.0,8.0) /), shape(v))
  call zvecprint(n, v, "v")
  call zgeev('N', '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(n, v, "lambda")
  call zmatprint(n, m, vR, "values")
  call zcheckeigen(n, A, e, vR);

end program test_eigen

! 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

  print *, "zCheckEigen"
  do k = 1,n
    do j = 1,n
      vk(j) = V(k,j)
      ve(j) = e(k)*V(k,j)
    end do
    call zMulMatVec(n, n, A, vk, Av);
    call zvecprint(n, vk, "v ")
    call zvecprint(n, ve, "ev")
    call zvecprint(n, Av, "Av")
    do m = 1,n
      Avve(m) = abs(Av(m)-ve(m))
    end do
    call zvecprint(n, 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(n, v, c)
  implicit none
  integer, intent(in) :: n
  double complex, dimension(n), intent(in) :: v
  character, intent(in) :: c
  integer :: i
  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.
!