Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/examples/optim1.f90
03540012fedefeature/fft-split

Last Change
Raw File

optim1.f90
View Options

!>
!> @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/>.
!>
!> @authors
!> (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
c4science · Help