!
!
!
program unit_test0

  use dc_types, only: &
    DP

  use dc_test, only: AssertLessThan, AssertEqual

  use igmcore_math, only: &
    & PI

  use igmcore_SimLinearAlgebraEq, only: &
    & LUdecomp_TriDiagMat, sol_SimLinearAlEq_lutri

  use igmcore_Interpolation, only: &
    & CubicSpline_Interpolation, CubicSpline_eval, &
    & CubicSpline_SData

  implicit none

  ! Work variables for unit test 0.1 and 0.2
  !
  integer, parameter :: N = 4
  integer :: i,j

  real(DP) :: al(N), ad(N), au(N), b(N)
  real(DP) :: error, exanswervec(N), answerVec(N)

  ! WOrk variables for unit test 0.3
  !
  integer, parameter :: CSDATA_N = 50
  real(DP) :: csGrid(CSDATA_N), csdata(CSDATA_N)
  type(CubicSpline_SData) :: csS(CSDATA_N-1)
  real(DP) :: check_point, check_val

  ! Unit test 0.1
  !
  write(*,*) '= Unit Test 0.1 [ igmcore_SimLinearAlgebraEq ]'
  write(*,*) '* Check the solver of simultaneous linear algebraic equation containg tridiagonal matrix.'

  write(*,*) '-- construct tridiagonal matrix ..'

  al(:) = -1.0d0
  ad(:) = 2.0d0
  au(:) = -1.0d0
  b(1:N) = (/ 1.0d0, 0.0d0, 0.0d0, 1.0d0 /)

  do i=1, N
    write(*,*) al(i), ad(i), au(i)
  end do

  write(*,*) '-- construct RHS vector ..'
  write(*,*) b(1:N)

  call LUdecomp_TriDiagMat(al(:), ad(:), au(:))
  write(*,*) '-- LU decomposition of tridiagonal matrix ..'
  do i=1, N
    write(*,*) al(i), ad(i), au(i)
  end do


  write(*,*) '-- solve simultaneous linear algebraic equation ..'

  answerVec(:) = sol_SimLinearAlEq_lutri(al(:), ad(:), au(:), b(:))
  write(*,*) 'answer vector=', answerVec(1:N)

  exanswervec(:) = 1.0d0
  error = dot_product(answerVec(:)-exanswervec(:), answerVec(:)-exanswervec(:))

  call AssertLessThan( message='Check answer', &
    & answer=1.0d-12, check=error  )


  ! Unit test 0.2
  !
  write(*,*) '= Unit Test 0.2 [ igmcore_SimLinearAlgebraEq ]'
  write(*,*) '* Check the solver of simultaneous linear algebraic equation containg tridiagonal matrix.'

  write(*,*) '-- construct tridiagonal matrix ..'

  al(:) = -1.0d0
  ad(:) = 2.0d0
  ad(N) = 1.0d0
  au(:) = -1.0d0
  b(1:N) = (/ 0.0d0, 0.0d0, 0.0d0, 1.0d0 /)

  do i=1, N
    write(*,*) al(i), ad(i), au(i)
  end do

  write(*,*) '-- construct RHS vector ..'
  write(*,*) b(1:N)

  call LUdecomp_TriDiagMat(al(:), ad(:), au(:))
  write(*,*) '-- LU decomposition of tridiagonal matrix ..'
  do i=1, N
    write(*,*) al(i), ad(i), au(i)
  end do


  write(*,*) '-- solve simultaneous linear algebraic equation ..'

  answerVec(:) = sol_SimLinearAlEq_lutri(al(:), ad(:), au(:), b(:))
  write(*,*) 'answer vector=', answerVec(1:N)
	
	do i=1,N
  	exanswervec(i) = i
  end do
  
  error = dot_product(answerVec(:)-exanswervec(:), answerVec(:)-exanswervec(:))

  call AssertLessThan( message='Check answer', &
    & answer=1.0d-12, check=error  )

  ! Unit test 0.3
  !
  write(*,*) '= Unit Test 0.3 [ igmcore_Interpolation ]'
  write(*,*) '* Check a procedure for interpolation using the cubic spline interpolation.'

  do i=1, CSDATA_N
    csGrid(i) = 2d0 * PI / real( CSDATA_N - 1 ) * ( i - 1 )
    csdata(i) = sin(csGrid(i))
  end do

  csS(:) = CubicSpline_Interpolation(csGrid, csdata)

  write(*,*) '-- pos, calculated value, answer'

  error = 0d0
  do i=1, CSDATA_N-1
    check_point = csGrid(i) + PI / real( CSDATA_N - 1 )
    check_val = CubicSpline_eval(csS(i), check_point)
    write(*,*) check_point, check_val, sin(check_point)
    error = error + abs(check_val - sin(check_point)) / sin(check_point)
  end do

  error = error / CSDATA_N * 100d0
  write(*,*) 'relative error(%): ', error

  ! Check whether relative error for the interpolaton is less than 0.1 %.
  call AssertLessThan( message='Check accuracy for cubic spline interpolation', &
    & answer=1.0d-1, check=error )

end program unit_test0
