!>
!> @file optim1.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   Test and compare performance of using "spline" and
!   "pp" forms. 1D case
!
  USE bsplines
!
  IMPLICIT NONE
  INTEGER :: nx, nidbas, nrank, npt=1000000
  INTEGER :: i
  DOUBLE PRECISION :: dx
  DOUBLE PRECISION :: seconds, t0, t1
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, coefs
  DOUBLE PRECISION, ALLOCATABLE :: xpt(:), fpt0(:), fpt1(:), fun(:, :)
  INTEGER, ALLOCATABLE :: left(:)
  TYPE(spline1d) :: splx
!
  NAMELIST /newrun/ nx, nidbas, npt
!===============================================================================
!
!   1D grid
!
  nx = 10
  nidbas = 3
  npt = 1000000
  READ(*,newrun)
  WRITE(*,newrun)

  ALLOCATE(xgrid(0:nx))
  dx = 1.0d0/REAL(nx)
  xgrid = (/ (i*dx,i=0,nx) /)
  WRITE(*,'(a/(10f8.3))') 'xgrid', xgrid
!
!   Set up spline
!
  CALL set_spline(nidbas, 4, xgrid, splx)
  nrank = splx%dim
  WRITE(*,'(a, i5)') 'nrank =', nrank
  WRITE(*,'(a/(10f8.3))') 'knots', splx%knots
!
  ALLOCATE(xpt(npt))
  ALLOCATE(left(npt))
  ALLOCATE(fun(0:nidbas,0:1))  ! Values and first derivatives of all Splines
  CALL RANDOM_NUMBER(xpt)
  CALL locintv(splx, xpt, left)
!===============================================================================
!
!   Check def_basfun_opt
!
  CALL basfun_recur(xpt(101), splx, fun, left(101)+1)
  WRITE(*,'(/a,f20.15,i4/(2f20.15))') 'BASFUN_RECUR at X=', xpt(101), left(101),&
       & (fun(:,i), i=0,1)
!
!!$  CALL def_basfun(xpt(101), splx, fun)
  CALL basfun(xpt(101), splx, fun, left(101)+1)
  WRITE(*,'(/a,f20.15/(2f20.15))') 'DEF_BASFUN at X=', xpt(101), &
       & (fun(:,i), i=0,1)
!
!   Performance of basis function computations
!
  t0 = seconds()
  DO i=1,npt
     CALL basfun_recur(xpt(i), splx, fun, left(i)+1)
  END DO
  WRITE(*,'(/a,1pe12.3)') 'BASFUN_RECUR time (s)', (seconds()-t0)/REAL(npt)
!
  t0 = seconds()
  DO i=1,npt
!!$     CALL def_basfun(xpt(i), splx, fun)
     CALL basfun(xpt(i), splx, fun, left(i)+1)
  END DO
  WRITE(*,'(/a,1pe12.3)') 'DEF_BASFUN time (s)', (seconds()-t0)/REAL(npt)
!===============================================================================
!
!   Check and performance of GRIDVAL using PPFORM and SPLINE expansion
!
  ALLOCATE(coefs(nrank))
  DEALLOCATE(xpt)
  ALLOCATE(xpt(npt), fpt0(npt), fpt1(npt))
  CALL RANDOM_NUMBER(xpt)
!
  splx%nlppform = .TRUE.
  coefs = 1.0d0
!
  CALL gridval(splx, xpt(1:1), fpt0(1:1), 0, coefs)
!
  t0 = seconds()
  CALL gridval(splx, xpt, fpt1, 1)
  CALL gridval(splx, xpt, fpt0, 0)
  t1 = seconds()-t0
  WRITE(*,'(/a,2(1pe12.3))') 'PPFORM: Max errors', &
       & MAXVAL(ABS(fpt0-1.0d0)) ,MAXVAL(ABS(fpt1))
  WRITE(*,'(a,3(1pe12.3))') 'time(s)', t1/REAL(npt)
!
  splx%nlppform = .FALSE.
  coefs = 1.0d0
  CALL gridval(splx, xpt(1:1), fpt0(1:1), 0, coefs)
  t0 = seconds()
  CALL gridval(splx, xpt, fpt1, 1)
  CALL gridval(splx, xpt, fpt0, 0)
  t1 = seconds()-t0
  WRITE(*,'(/a,2(1pe12.3))') 'SPLINES: Max errors', &
       & MAXVAL(ABS(fpt0-1.0d0)) ,MAXVAL(ABS(fpt1))
  WRITE(*,'(a,3(1pe12.3))') 'time(s)', t1/REAL(npt)
!===============================================================================
!
!   Clean up
!
  CALL destroy_sp(splx)
  DEALLOCATE(xgrid, coefs)
  DEALLOCATE(xpt, fpt0, fpt1)
  DEALLOCATE(fun)
END PROGRAM main