!>
!> @file optim3.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. 2D1D case
!
  USE bsplines
!
  IMPLICIT NONE
  INTEGER :: nx, ny, nz, ngauss(3), nidbas(3)
  INTEGER :: npt, d1, d2, d3
  INTEGER :: i, j, k
  DOUBLE PRECISION :: pi, dx, dy, dz
  DOUBLE PRECISION :: seconds, t0, t1
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, ygrid, zgrid
  DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: coefs
  DOUBLE PRECISION, ALLOCATABLE :: xpt(:), ypt(:), zpt(:)
  DOUBLE PRECISION, ALLOCATABLE :: fpt000(:), fpt100(:), fpt010(:), fpt001(:)
  TYPE(spline2d1d) :: splxyz
!
  NAMELIST /newrun/ nx, ny, nz, nidbas, npt
!===============================================================================
!
!   2D grid
!
  nx = 8
  ny = 8
  nz = 8
  nidbas = (/ 3, 3, 3 /)
  npt = 1000000
  READ(*,newrun)
  WRITE(*,newrun)
!
  pi = 4.0d0*ATAN(1.0d0)
  ALLOCATE(xgrid(0:nx), ygrid(0:ny), zgrid(0:nz))
  dx = 1.0d0/REAL(nx)
  xgrid = (/ (i*dx,i=0,nx) /)
  dy = 2.0d0*pi/REAL(ny)
  ygrid = (/ (j*dy,j=0,ny) /)
  dz = 2.0d0*pi/REAL(nz)
  zgrid = (/ (k*dz,k=0,nz) /)
!
  WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
  WRITE(*,'(a/(10f8.3))') 'YGRID', ygrid(0:ny)
  WRITE(*,'(a/(10f8.3))') 'ZGRID', zgrid(0:nz)
!
!   Set up spline
!
  ngauss = 4
  CALL set_spline(nidbas, ngauss, xgrid, ygrid, zgrid, splxyz, &
       &         (/.FALSE., .TRUE., .TRUE./))
  d1 = splxyz%sp12%sp1%dim
  d2 = splxyz%sp12%sp2%dim
  d3 = splxyz%sp3%dim
  WRITE(*,'(a,3i4)') 'd1, d2, d3 =', d1, d2, d3
  WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in X', splxyz%sp12%sp1%knots
  WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in Y', splxyz%sp12%sp2%knots
  WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in Z', splxyz%sp3%knots
!===============================================================================
!
!   Check and performance of GRIDVAL using PPFORM and SPLINE expansion
!
  ALLOCATE(xpt(npt), ypt(npt), zpt(npt))
  CALL RANDOM_NUMBER(xpt)
  CALL RANDOM_NUMBER(ypt); ypt = 2.0d0*pi*ypt
  CALL RANDOM_NUMBER(zpt); zpt = 2.0d0*pi*zpt
!
  ALLOCATE(coefs(d1,d2,d3))
  ALLOCATE(fpt000(npt), fpt100(npt), fpt010(npt), fpt001(npt))
!
  coefs = 1.0d0
!
  splxyz%sp12%sp1%nlppform = .TRUE.
  splxyz%sp12%sp2%nlppform = .TRUE.
  splxyz%sp3%nlppform = .TRUE.
  CALL gridval(splxyz, xpt, ypt, zpt, fpt000, (/0,0,0/), coefs)
  t0 = seconds()
  CALL gridval(splxyz, xpt, ypt, zpt, fpt000, (/0,0,0/))
  CALL gridval(splxyz, xpt, ypt, zpt, fpt100, (/1,0,0/))
  CALL gridval(splxyz, xpt, ypt, zpt, fpt010, (/0,1,0/))
  CALL gridval(splxyz, xpt, ypt, zpt, fpt001, (/0,0,1/))
  t1 = seconds()-t0
  WRITE(*,'(/a,4(1pe12.3))') 'PPFORM: Max errors', &
       & MAXVAL(ABS(fpt000-1.0d0)), &
       & MAXVAL(ABS(fpt100)), &
       & MAXVAL(ABS(fpt010)), &
       & MAXVAL(ABS(fpt001))
  WRITE(*,'(a,3(1pe12.3))') 'time(s)', t1/REAL(npt)
!
  splxyz%sp12%sp1%nlppform = .FALSE.
  splxyz%sp12%sp2%nlppform = .FALSE.
  splxyz%sp3%nlppform = .FALSE.
  CALL gridval(splxyz, xpt, ypt, zpt, fpt000, (/0,0,0/), coefs)
  t0 = seconds()
  CALL gridval(splxyz, xpt, ypt, zpt, fpt000, (/0,0,0/))
  CALL gridval(splxyz, xpt, ypt, zpt, fpt100, (/1,0,0/))
  CALL gridval(splxyz, xpt, ypt, zpt, fpt010, (/0,1,0/))
  CALL gridval(splxyz, xpt, ypt, zpt, fpt001, (/0,0,1/))
  t1 = seconds()-t0
  WRITE(*,'(/a,4(1pe12.3))') 'SPLINES: Max errors', &
       & MAXVAL(ABS(fpt000-1.0d0)), &
       & MAXVAL(ABS(fpt100)), &
       & MAXVAL(ABS(fpt010)), &
       & MAXVAL(ABS(fpt001))
  WRITE(*,'(a,3(1pe12.3))') 'time(s)', t1/REAL(npt)
!===============================================================================
!
!   Clean up
!
  CALL destroy_sp(splxyz)
  DEALLOCATE(xgrid, ygrid, zgrid, coefs)
  DEALLOCATE(xpt, ypt, fpt000,fpt100, fpt010, fpt001)
END PROGRAM main