!>
!> @file tlocintv.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
!
!   Optimization of locintv
!
  USE bsplines
  IMPLICIT NONE
  INTEGER :: nx, nidbas, ngauss, np, nits
  DOUBLE PRECISION :: a, b, coefs(5)
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid
  TYPE(spline1d) :: splx
!
  INTEGER :: i, nerrs, it
  DOUBLE PRECISION :: t0, t1, seconds, tscal, tscal_new, tvec, tvec_new
  DOUBLE PRECISION, ALLOCATABLE :: xp(:)
  INTEGER, ALLOCATABLE :: left(:)
!
  INTERFACE
     SUBROUTINE meshdist(coefs, x, nx)
       DOUBLE PRECISION, INTENT(in) :: coefs(5)
       INTEGER, INTENT(iN) :: nx
       DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
     END SUBROUTINE meshdist
  END INTERFACE
!
  NAMELIST /newrun/ nx, nidbas, ngauss, np, nits, a, b, coefs
!===========================================================================
!   Read in data
!
  nx = 8       ! Number oh intevals in x
  nidbas = 3   ! Degree of splines
  ngauss = 4   ! Number of Gauss points/interval
  np = 10 ! Number of random points in [a,b]
  nits = 1000000
  a = 0.0
  b = 1.0
  coefs(1:5) = (/0.2d0, 1.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
                                                 ! see function FDIST in MESHDIST
!
  READ(*,newrun)
  WRITE(*,newrun)
!
!   Define grid on x axis
!
  ALLOCATE(xgrid(0:nx))
  xgrid(0) = a
  xgrid(nx) = b
  CALL meshdist(coefs, xgrid, nx)
  WRITE(*,'(/a/(10f8.3))') 'XGRID', xgrid(0:nx)
  WRITE(*,'(a/(10f8.3))') 'Diff XGRID', (xgrid(i)-xgrid(i-1), i=1,nx)
!
!   Set up spline
!
  CALL set_spline(nidbas, ngauss, xgrid, splx)
  WRITE(*,'(a,l1)') 'Is mesh equidistant? ', splx%nlequid
!
!   Test locintv
!
  ALLOCATE(xp(np))
  ALLOCATE(left(np))
  xp(:) = (b-a)*xp(:) + a
  tscal = 0.0d0
  tscal_new = 0.0d0
  tvec = 0.0d0
  tvec_new = 0.0d0
!
  nerrs = 0
  DO it=1,nits
     CALL RANDOM_NUMBER(xp)
     t0 = seconds()
     DO i=1,np
        CALL locintv_old(splx, xp(i), left(i))
     END DO
     tscal = tscal + seconds()-t0
     nerrs = nerrs + COUNT(.NOT.in_interv(xp, left))
!
     t0 = seconds()
     DO i=1,np
        CALL locintv(splx, xp(i), left(i))
     END DO
     tscal_new = tscal_new + seconds()-t0
     nerrs = nerrs + COUNT(.NOT.in_interv(xp, left))
!
     t0 = seconds()
     CALL locintv_old(splx, xp, left)
     tvec = tvec + seconds()-t0
     nerrs = nerrs + COUNT(.NOT.in_interv(xp, left))
!
     t0 = seconds()
     CALL locintv(splx, xp, left)
     tvec_new = tvec_new + seconds()-t0
     nerrs = nerrs + COUNT(.NOT.in_interv(xp, left))
  END DO
  PRINT*, 'nerrs =', nerrs
!
  tscal = tscal/(REAL(nits*np,8))
  tscal_new = tscal_new/(REAL(nits*np,8))
  tvec = tvec/(REAL(nits*np,8))
  tvec_new = tvec_new/(REAL(nits*np,8))
  WRITE(*,'(4a12)') 'scalar', 'scalar new', 'vector', 'vector new'
  WRITE(*,'(4(1pe12.3))') tscal, tscal_new, tvec, tvec_new
!
!   Clean up
!
  DEALLOCATE(xp)
  DEALLOCATE(left)
  DEALLOCATE(xgrid)
  CALL destroy_sp(splx)
!
CONTAINS
  LOGICAL ELEMENTAL FUNCTION in_interv(x, l)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(in) :: l
    in_interv = x.GE.xgrid(l) .AND. x.LT.xgrid(l+1)
  END FUNCTION in_interv  
END PROGRAM main