Homec4science
Diffusion SPClibs (master) Browse

spclibs/examples/gridval_perf.f90
c723281d03b6master

Last Change
Raw File

gridval_perf.f90
View Options

!>
!> @file gridval_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/>.
!>
!> @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. 2D case
!
USE bsplines
!
IMPLICIT NONE
INTEGER :: nx, ny, ngauss(2), nidbas(2), nits
INTEGER :: npt, d1, d2
INTEGER :: i, j, its, ngroup=4
INTEGER :: i1, i2, nset, nremain
DOUBLE PRECISION :: pi, dx, dy
DOUBLE PRECISION :: seconds, t0, t1, tscal, tvec
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, ygrid
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: coefs
DOUBLE PRECISION, ALLOCATABLE :: xpt(:), ypt(:), fpt00(:), fpt01(:), fpt10(:)
DOUBLE PRECISION, ALLOCATABLE :: fscal00(:), fscal01(:), fscal10(:)
TYPE(spline2d) :: splxy
!
NAMELIST /newrun/ nx, ny, nidbas, npt, nits
!===============================================================================
! 0.0 Prologue
!
! 2D grid
!
nx = 8
ny = 8
nidbas = (/ 3, 3 /)
npt = 100000
nits =100
READ(*,newrun)
WRITE(*,newrun)
!
pi = 4.0d0*ATAN(1.0d0)
ALLOCATE(xgrid(0:nx), ygrid(0:ny))
dx = 1.0d0/REAL(nx)
xgrid = (/ (i*dx,i=0,nx) /)
dy = 2.0d0*pi/REAL(ny)
ygrid = (/ (j*dy,j=0,ny) /)
!
WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
WRITE(*,'(a/(10f8.3))') 'YGRID', ygrid(0:ny)
!
! Set up spline
!
ngauss = 4
CALL set_spline(nidbas, ngauss, xgrid, ygrid, splxy, (/.FALSE., .TRUE./))
d1 = splxy%sp1%dim
d2 = splxy%sp2%dim
WRITE(*,'(a,3i4)') 'd1, d2 =', d1, d2
WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in X', splxy%sp1%knots
WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in Y', splxy%sp2%knots
!===============================================================================
! 1.0 PPFORM
!
ALLOCATE(xpt(npt), ypt(npt))
CALL RANDOM_NUMBER(xpt)
CALL RANDOM_NUMBER(ypt); ypt = 2.0d0*pi*ypt
!
ALLOCATE(coefs(d1,d2))
ALLOCATE(fscal00(npt), fscal01(npt), fscal10(npt))
ALLOCATE(fpt00(npt), fpt01(npt), fpt10(npt))
!
coefs = 1.0d0 ! => f=1, all derivatives = 0!
!
splxy%sp1%nlppform = .TRUE.
splxy%sp2%nlppform = .TRUE.
CALL gridval(splxy, xpt, ypt, fscal00, (/0,0/), coefs)
!
! Scalar PPFORM
t0 = seconds()
DO i=1,npt
DO its=1,nits
CALL gridval(splxy, xpt(i), ypt(i), fscal00(i), (/0,0/))
CALL gridval(splxy, xpt(i), ypt(i), fscal01(i), (/0,1/))
CALL gridval(splxy, xpt(i), ypt(i), fscal10(i), (/1,0/))
END DO
END DO
t1 = seconds()-t0
tscal = t1/REAL(npt*nits,8)
WRITE(*,'(/a,3(1pe12.3))') 'Scalar PPFORM: Max errors', &
& MAXVAL(ABS(fscal00-1.0d0)), MAXVAL(ABS(fscal01)), MAXVAL(ABS(fscal10))
WRITE(*,'(a,3(1pe12.3))') 'time(s)', tscal
!
! Vector PPFORM
WRITE(*,'(/a/a5,2a12)') 'Vector ppform', 'N', 't(s)', 'SpeedUp'
ngroup = 1
DO WHILE (ngroup.LT.npt/2)
nset = npt/ngroup
nremain = MODULO(npt, ngroup)
IF(nremain.NE.0) nset = nset+1
t0 = seconds()
i2 = 0
DO i=1,nset
i1 = i2+1
i2 = MIN(i2+ngroup,npt)
DO its=1,nits
CALL gridval(splxy, xpt(i1:i2), ypt(i1:i2), fpt00(i1:i2), (/0,0/))
CALL gridval(splxy, xpt(i1:i2), ypt(i1:i2), fpt01(i1:i2), (/0,1/))
CALL gridval(splxy, xpt(i1:i2), ypt(i1:i2), fpt10(i1:i2), (/1,0/))
END DO
END DO
t1 = seconds()-t0
tvec = t1/REAL(npt*nits)
WRITE(*,'(i5,3(1pe12.3))') ngroup, tvec, tscal/tvec
ngroup = 2*ngroup
END DO
WRITE(*,'(/a,3(1pe12.3))') 'Vector PPFORM: Max errors', &
& MAXVAL(ABS(fpt00-fscal00)), MAXVAL(ABS(fpt01-fscal01)), MAXVAL(ABS(fpt10-fscal10))
!===============================================================================
! 2.0 Sline expansion
!
coefs = 1.0d0 ! => f=1, all derivatives = 0!
!
splxy%sp1%nlppform = .FALSE.
splxy%sp2%nlppform = .FALSE.
CALL gridval(splxy, xpt, ypt, fscal00, (/0,0/), coefs)
!
! Scalar SPLINE
t0 = seconds()
DO i=1,npt
DO its=1,nits
CALL gridval(splxy, xpt(i), ypt(i), fscal00(i), (/0,0/))
CALL gridval(splxy, xpt(i), ypt(i), fscal01(i), (/0,1/))
CALL gridval(splxy, xpt(i), ypt(i), fscal10(i), (/1,0/))
END DO
END DO
t1 = seconds()-t0
tscal = t1/REAL(npt*nits,8)
WRITE(*,'(/a,3(1pe12.3))') 'Scalar SPLINE: Max errors', &
& MAXVAL(ABS(fscal00-1.0d0)), MAXVAL(ABS(fscal01)), MAXVAL(ABS(fscal10))
WRITE(*,'(a,3(1pe12.3))') 'time(s)', tscal
!
! Vector SPLINE
WRITE(*,'(/a/a5,2a12)') 'Vector spline', 'N', 't(s)', 'SpeedUp'
ngroup = 1
DO WHILE (ngroup.LT.npt/2)
nset = npt/ngroup
nremain = MODULO(npt, ngroup)
IF(nremain.NE.0) nset = nset+1
t0 = seconds()
i2 = 0
DO i=1,nset
i1 = i2+1
i2 = MIN(i2+ngroup,npt)
DO its=1,nits
CALL gridval(splxy, xpt(i1:i2), ypt(i1:i2), fpt00(i1:i2), (/0,0/))
CALL gridval(splxy, xpt(i1:i2), ypt(i1:i2), fpt01(i1:i2), (/0,1/))
CALL gridval(splxy, xpt(i1:i2), ypt(i1:i2), fpt10(i1:i2), (/1,0/))
END DO
END DO
t1 = seconds()-t0
tvec = t1/REAL(npt*nits)
WRITE(*,'(i5,3(1pe12.3))') ngroup, tvec, tscal/tvec
ngroup = 2*ngroup
END DO
WRITE(*,'(/a,3(1pe12.3))') 'Vector SPLINE: Max errors', &
& MAXVAL(ABS(fpt00-fscal00)), MAXVAL(ABS(fpt01-fscal01)), MAXVAL(ABS(fpt10-fscal10))
!===============================================================================
!
! Clean up
!
CALL destroy_sp(splxy)
DEALLOCATE(xgrid, ygrid, coefs)
DEALLOCATE(xpt, ypt, fpt00, fpt01, fpt10)
DEALLOCATE(fscal00, fscal01, fscal10)
END PROGRAM main
c4science · Help