!------------------------------------------------------------------------
! Copyright (c) 2005-2011 SPMODEL Development Group. All rights reserved.
!------------------------------------------------------------------------
!
!表題  ceigen/eigmatrix テストプログラム
!
!履歴  2005/01/26  竹広真一
!      2007/11/02  竹広真一  メッセージ, 正誤判定追加
!      2008/12/05  佐々木洋平. 文字列 sort の長さを 2 に
!      2011/02/18  佐々木洋平 dc_test を使用するように修正
!      2011/11/05  竹広真一  ceigen 用に改造
!      2019/10/10  佐々木洋平 module rename
!
program ceigen_test

  use dc_message, only : MessageNotify
  use dc_types, only: DP
  use dc_test, only : AssertEqual
  use spml_utils, only: ceigen
  implicit none

  complex(DP), dimension(:,:), allocatable  :: cmatrix
  complex(DP), dimension(:),   allocatable  :: eigval
  complex(DP), dimension(:,:), allocatable  :: eigvec

  complex(DP), dimension(:),   allocatable       :: eigval_sol
  complex(DP), dimension(:,:), allocatable       :: eigvec_sol
  integer                     :: info
  integer                     :: i
  complex(DP), parameter      :: EI = (0.0D0, 1.0D0)
  ! 判定誤差設定
  integer, parameter :: check_digits = 10
  integer, parameter :: ignore = -11

  call MessageNotify('M','ceigen_test','Test of ceigen/eigmatrix')

!--------------- 2x2 行列 ------------------
!  VALUE = (0,1), (0,-1)
!  VECTOR = (-1,1), (1,1)

  allocate(cmatrix(2,2))
  allocate(eigval(2),eigval_sol(2))
  allocate(eigvec(2,2),eigvec_sol(2,2))


  !---- 行列 ----
  cmatrix(1,:) = (/ (0.0d0,0.0D0),     EI/)
  cmatrix(2,:) = (/ EI,     (0.0d0,0.0D0)/)

  !---- 正解 ----
  eigval_sol = (/EI,-EI/)

  eigvec_sol(:,1) = (/(1.0d0,0.0d0),(1.0d0,0.0d0)/)
  eigvec_sol(:,2) = (/(-1.0d0,0.0d0),(1.0d0,0.0d0)/)

  !---- 固有値計算 ----
  call ceigen(cmatrix,eigval,eigvec,info,sort='I ',reverse=.true.)
  do i=1,2
     eigvec(:,i) = eigvec(:,i)/eigvec(1,i)
     eigvec_sol(:,i) = eigvec_sol(:,i)/eigvec_sol(1,i)
  enddo

  !---- 判定 ----
  call AssertEqual(                                           &
    message= '2x2 arrray eigenvalue(real part)',              &
    answer = dble(eigval_sol),                                &
    check = dble(eigval),                                     &
    significant_digits = check_digits, ignore_digits = ignore &
    )

  call AssertEqual(                                           &
    message= '2x2 arrray eigenvalue(imaginary part)',         &
    answer = aimag(eigval_sol),                               &
    check = aimag(eigval),                                    &
    significant_digits = check_digits, ignore_digits = ignore &
    )

  call AssertEqual(                                           &
    message= '2x2 arrray eigenvector(real part)',             &
    answer = dble(eigvec_sol),                                &
    check = dble(eigvec),                                     &
    significant_digits = check_digits, ignore_digits = ignore &
    )

  call AssertEqual(                                           &
    message= '2x2 arrray eigenvector(imaginary part)',        &
    answer = aimag(eigvec_sol),                               &
    check = aimag(eigvec),                                    &
    significant_digits = check_digits, ignore_digits = ignore &
    )

  call MessageNotify('M','Test of ceigen/eigmatrix','Test of 2x2 array succeeded!')

  deallocate(cmatrix)
  deallocate(eigval,eigval_sol)
  deallocate(eigvec,eigvec_sol)

!--------------- 3x3 行列 ------------------
!  VALUE = (6,0), (3,0), (2,0)
!  VECTOR = (-1,2i,1), (1,i,-1),(1,0,1)

  allocate(cmatrix(3,3))
  allocate(eigval(3), eigval_sol(3))
  allocate(eigvec(3,3), eigvec_sol(3,3))

  !---- 行列 ----
  cmatrix(1,:) = (/(3.0d0,0.0D0),       EI, (-1.0D0,0.0D0)/)
  cmatrix(2,:) = (/ -EI,     (5.0d0,0.0D0),             EI/)
  cmatrix(3,:) = (/(-1.0d0,0.0d0),     -EI,  (3.0d0,0.0D0)/)

  !---- 正解 ----
  eigval_sol = (/(6.0D0,0.0D0),(3.0d0,0.0d0),(2.0d0,0.0d0)/)

  eigvec_sol(:,1) = (/(-1.0d0,0.0d0),      2.0d0*EI,  (1.0d0,0.0d0)/)
  eigvec_sol(:,2) = (/ (1.0d0,0.0d0),            EI, (-1.0d0,0.0d0)/)
  eigvec_sol(:,3) = (/ (1.0d0,0.0d0), (0.0d0,0.0d0),  (1.0d0,0.0d0)/)

  !---- 固有値計算 ----
  call ceigen(cmatrix,eigval,eigvec,info,sort='R ',reverse=.true.)

  do i=1,3
     eigvec(:,i) = eigvec(:,i)/eigvec(1,i)
     eigvec_sol(:,i) = eigvec_sol(:,i)/eigvec_sol(1,i)
  enddo

  !---- 判定 ----
  call AssertEqual(                                           &
    message= '3x3 arrray eigenvalue(real part)',              &
    answer = dble(eigval_sol),                                &
    check = dble(eigval),                                     &
    significant_digits = check_digits, ignore_digits = ignore &
    )

  call AssertEqual(                                           &
    message= '3x3 arrray eigenvalue(imaginary part)',         &
    answer = aimag(eigval_sol),                               &
    check = aimag(eigval),                                    &
    significant_digits = check_digits, ignore_digits = ignore &
    )

  call AssertEqual(                                           &
    message= '3x3 arrray eigenvector(real part)',             &
    answer = dble(eigvec_sol),                                &
    check = dble(eigvec),                                     &
    significant_digits = check_digits, ignore_digits = ignore &
    )

  call AssertEqual(                                           &
    message= '3x3 arrray eigenvector(imaginary part)',        &
    answer = aimag(eigvec_sol),                               &
    check = aimag(eigvec),                                    &
    significant_digits = check_digits, ignore_digits = ignore &
    )

  call MessageNotify('M','Test of ceigen/eigmatrix','Test of 3x3 array succeeded!')

  deallocate(cmatrix)
  deallocate(eigval,eigval_sol)
  deallocate(eigvec,eigvec_sol)

end program ceigen_test
