Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/tlocintv.f90
03540012fedefeature/refactoring

Last Change
Raw File

tlocintv.f90
View Options

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