!------------------------------------------------------------------------
! 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 用に改造
!
!備考  lapack, blas ライブラリが必要. Debian/GNU Linux + Fujitsu frt ならば
!      lapack, lapack-deb パッケージをインストールして,
!         -llapack -lblas -L/usr/lib/gcc-lib/i386-linux/2.95.4 -lg2c
!      といったオプションをつけるべし.
!
!
program ceigen_test

  use dc_message, only: MessageNotify
  use dc_test, only: AssertEqual
  use eigmatrix
  implicit none

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

  complex(8), dimension(:), allocatable :: eigval_sol
  complex(8), dimension(:, :), allocatable :: eigvec_sol
  integer :: info
  integer :: i
  complex(8), 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)
  end do

  !---- 判定 ----
  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)
  end do

  !---- 判定 ----
  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