!>
!> @file basfun_perf.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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   Performance of scalar and vector versions of def_basfun
!
  USE bsplines
  IMPLICIT NONE
  INTEGER :: nx, nidbas, nrank, npt=10, jdermx
  DOUBLE PRECISION :: dx
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid
  DOUBLE PRECISION, ALLOCATABLE :: xpt(:), fun(:, :)
  INTEGER :: left, i, i1, i2
  INTEGER :: ngroup, nset, nremain
  TYPE(spline1d) :: splx
  DOUBLE PRECISION :: t0, t1, seconds
  DOUBLE PRECISION :: t_loop, t_locintv1, t_basfun1, t_locintv, t_basfun
  INTEGER :: its, nits
  INTEGER, ALLOCATABLE :: vleft(:)
  DOUBLE PRECISION, ALLOCATABLE :: vfun(:,:,:)
  LOGICAL :: nlperiod
!
  NAMELIST /newrun/ nx, nidbas, npt, nits, ngroup, jdermx, nlperiod
!
!===============================================================================
!
!   1D grid
!
  nx = 10
  nidbas = 3
  npt = 1000000
  nits = 100
  ngroup = 10
  jdermx = 0
  nlperiod = .FALSE.
  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, period=nlperiod)
  nrank = splx%dim
  WRITE(*,'(a, i5)') 'nrank =', nrank
  WRITE(*,'(a/(10f8.3))') 'knots', splx%knots
!
  ALLOCATE(xpt(npt))
  ALLOCATE(fun(0:nidbas,0:jdermx))  ! Values and first derivatives of all Splines
  CALL RANDOM_NUMBER(xpt)
!===============================================================================
!                            1.0 Scalar version
!
!    loop
  t0 = seconds()
  DO its=1,nits
     DO i=1,npt
     END DO
  END DO
  t_loop = (seconds()-t0)/REAL(nits*npt,8)
!
!   locintv
  t0 = seconds()
  DO its=1,nits
     DO i=1,npt
        CALL locintv(splx, xpt(i), left)
     END DO
  END DO
  t_locintv1 = (seconds()-t0)/REAL(nits*npt,8)
!
!   def_basfun
  t0 = seconds()
  DO its=1,nits
     DO i=1,npt
        CALL locintv(splx, xpt(i), left)
        CALL basfun(xpt(i), splx, fun, left+1)
     END DO
  END DO
  t_basfun1 = (seconds()-t0)/REAL(nits*npt,8)
!
  WRITE(*,'(6x,3a12)') 'loop', 'locintv', 'basfun'
  WRITE(*,'(6x,8(1pe12.3))') t_loop, t_locintv1, t_basfun1
!===============================================================================
!                            2.0 Vector version
!
  ngroup = 1
  DO WHILE (ngroup .LT. npt/2)
     ALLOCATE(vleft(ngroup))
     ALLOCATE(vfun(0:nidbas, 0:jdermx, ngroup))
     nset = npt/ngroup
     nremain = MODULO(npt, ngroup)
     IF(nremain.NE.0) nset = nset+1
!
!    loop
     t0 = seconds()
     DO its=1,nits
        i2 = 0
        DO i=1,nset
           i1 = i2+1
           i2 = MIN(i2+ngroup,npt)
        END DO
     END DO
     t_loop = (seconds()-t0)/REAL(nits*nset,8)
!
!   locintv
     t0 = seconds()
     DO its=1,nits
        i2 = 0
        DO i=1,nset
           i1 = i2+1
           i2 = MIN(i2+ngroup,npt)
           CALL locintv(splx, xpt(i1:i2), vleft)
        END DO
     END DO
     t_locintv = (seconds()-t0)/REAL(nits*npt,8)
!
!   basfun
     t0 = seconds()
     DO its=1,nits
        i2 = 0
        DO i=1,nset
           i1 = i2+1
           i2 = MIN(i2+ngroup,npt)
           CALL locintv(splx, xpt(i1:i2), vleft)
           CALL basfun(xpt(i1:i2), splx, vfun, vleft+1)
        END DO
     END DO
     t_basfun = (seconds()-t0)/REAL(nits*npt,8)
!
     WRITE(*,'(i6,8(1pe12.3))') ngroup, t_loop, t_locintv, t_basfun, &
          &      t_locintv1/t_locintv, t_basfun1/t_basfun
     DEALLOCATE(vleft)
     DEALLOCATE(vfun)
     ngroup = ngroup*2
  END DO
!===============================================================================
!
!   Clean up
!
  CALL destroy_sp(splx)
  DEALLOCATE(xgrid)
  DEALLOCATE(xpt)
  DEALLOCATE(fun)
END PROGRAM main