diff --git a/src/bsplines.f90 b/src/bsplines.f90
index 13bf402..16fcdf0 100644
--- a/src/bsplines.f90
+++ b/src/bsplines.f90
@@ -1,4285 +1,4285 @@
!>
!> @file bsplines.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 Stephan Brunner <stephan.brunner@epfl.ch>
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
MODULE bsplines
!
! BSPLINES: A module to construct B-Splines of any order on
! non-equidistant mesh. Can be used for interpolation and
! Finite Element discretization.
!
! T.M. Tran, S. Brunner, CRPP-EPFL
! February 2007
!
USE matrix
IMPLICIT NONE
PRIVATE
PUBLIC :: spline1d, set_spline, get_dim, get_gauss, gridval
PUBLIC :: spline2d, spline2d1d, def_knots, allsplines
PUBLIC :: set_splcoef, get_splcoef
PUBLIC :: fintg, calc_integ, destroy_sp
PUBLIC :: gauleg, CompMassMatrix
PUBLIC :: basfun_recur, basfun, def_basfun, is_equid, locintv_old, locintv
PUBLIC :: calc_fftmass, calc_fftmass_old
PUBLIC :: init_dft, ft_basfun
PUBLIC :: getgrad
PUBLIC :: dftmap
!
DOUBLE PRECISION, PARAMETER :: pi=3.1415926535897931d0
!
TYPE dftmap
INTEGER :: n ! Total number of modes
INTEGER :: kmin, kmax ! Define the Foyrier window
INTEGER :: dk ! Number of modes in window
DOUBLE PRECISION :: dx ! Interval in real space
INTEGER, POINTER :: mode_couplings(:) => NULL() ! Table of mode couplings
DOUBLE COMPLEX, POINTER :: coefs(:,:) => NULL() ! The restricted Fourier coefs
END TYPE dftmap
!
TYPE spline1d
INTEGER :: order ! Spline order = spline degree + 1
INTEGER :: nints ! Number of knot intervals
INTEGER :: nsites ! Number of interpolation sites
INTEGER :: dim ! Dimension of spline space
INTEGER :: left=0 ! Save value used by routine LOCATE
LOGICAL :: period ! is the grid periodic ?
LOGICAL :: nlequid ! is the grid equidistant ?
LOGICAL :: nlppform ! Construct and build PPFORM in GRIDVAL if .TRUE.
DOUBLE PRECISION :: lperiod ! Periodicity
DOUBLE PRECISION :: hinv ! Inverse of mesh size, when nlequid=T
INTEGER, POINTER ::&
& fmap(:) => NULL() ! Mapping of fine to coarse mesh
DOUBLE PRECISION, POINTER :: &
& knots(:) => NULL(), & ! Spline knots (-p:n+p)
& val0(:,:,:) => NULL(), & ! Values and deriv. at left boundary of all splines
& valc(:,:) => NULL(), & ! Values and deriv. at left boundary of equid. per. splines
& gausx(:,:) => NULL(), & ! Gauss abscissas
& gausw(:,:) => NULL(), & ! Gauss weights
& intspl(:) => NULL(), & ! Integral of splines
& ppform(:,:) => NULL(), & ! PPFORM coefs
& bcoefs(:) => NULL() ! Spline coefs
DOUBLE COMPLEX, POINTER :: &
& ppformz(:,:) => NULL(), & ! PPFORM coefs for complex function
& bcoefsc(:) => NULL() ! Spline coefs for complex function
TYPE(GBMAT) :: mat ! Interpolation matrix
TYPE(periodic_mat) :: matp ! Interpolation matrix (periodic case)
TYPE(dftmap) :: dft ! Define DFT mapping
END TYPE spline1d
!
TYPE spline2d
TYPE(spline1d) :: sp1 ! Spline in direction 1
TYPE(spline1d) :: sp2 ! Spline in direction 2
DOUBLE PRECISION, POINTER :: ppform(:,:,:,:) => NULL() ! 2d PPFORM coefs
DOUBLE PRECISION, POINTER :: bcoefs(:,:) => NULL() ! Spline coefs
DOUBLE COMPLEX, POINTER :: ppformz(:,:,:,:) => NULL() ! PPFORM coefs for complex function
DOUBLE COMPLEX, POINTER :: bcoefsc(:,:) => NULL() ! Spline coefs for complex function
END TYPE spline2d
!
TYPE spline2d1d
TYPE(spline2d) :: sp12 ! 2D spline for dir. 1 and 2
TYPE(spline1d) :: sp3 ! 1D spline for dir. 3
DOUBLE PRECISION, POINTER :: ppform(:,:,:,:,:,:) => NULL() ! PPFORM coefs
DOUBLE PRECISION, POINTER :: bcoefs(:,:,:) => NULL() ! Spline coefs
DOUBLE PRECISION, POINTER :: ppformz(:,:,:,:,:,:) => NULL()! PPFORM coefs for complex function
DOUBLE COMPLEX, POINTER :: bcoefsc(:,:,:) => NULL() ! Spline coefs for complex function
END TYPE spline2d1d
!
INTERFACE set_spline
MODULE PROCEDURE set_spline1d, set_spline2d, set_spline2d1d
END INTERFACE
INTERFACE get_dim
MODULE PROCEDURE get_dim1, get_dim2
END INTERFACE
INTERFACE gridval
MODULE PROCEDURE gridval1d, gridval1dz, &
& gridval2d, gridval2dz, &
& gridval2d_1d, gridval2d_1dz, &
& gridval2d_2d, gridval2d_2dz, &
& gridval2d1d_3d, gridval2d1d_1d
END INTERFACE
INTERFACE set_splcoef
MODULE PROCEDURE set_splcoef1d, set_splcoef2d
END INTERFACE
INTERFACE get_splcoef
MODULE PROCEDURE get_splcoef1, get_splcoef1z, get_splcoefn, &
& get_splcoef2d, get_splcoef2dz
END INTERFACE
INTERFACE fintg
MODULE PROCEDURE fintg1, fintg2
END INTERFACE
INTERFACE destroy_sp
MODULE PROCEDURE destroy_sp1d, destroy_sp2d, destroy_sp2d1d
END INTERFACE
INTERFACE CompMassMatrix
MODULE PROCEDURE CompMassMatrix1, CompMassMatrix_gb, CompMassMatrix_zgb
END INTERFACE
INTERFACE calc_integ
MODULE PROCEDURE calc_integ0,calc_integn
END INTERFACE
INTERFACE locintv
MODULE PROCEDURE locintv0, locintv1
END INTERFACE locintv
INTERFACE locintv_old
MODULE PROCEDURE locintv0_old, locintv1_old
END INTERFACE locintv_old
INTERFACE ppval
MODULE PROCEDURE ppval0, ppval1, ppval2, &
& ppval0_n, &
& ppval0z, ppval1z, ppval2z, &
& ppval0z_n
END INTERFACE ppval
INTERFACE basfun
MODULE PROCEDURE basfun0, basfun1
END INTERFACE basfun
INTERFACE ft_basfun
MODULE PROCEDURE ft_basfun0, ft_basfun1
END INTERFACE ft_basfun
INTERFACE def_basfun
MODULE PROCEDURE def_basfun0, def_basfun1
END INTERFACE def_basfun
INTERFACE getgrad
MODULE PROCEDURE getgradr, getgradz
END INTERFACE getgrad
!
CONTAINS
!===========================================================================
SUBROUTINE set_spline1d(p, ngauss, grid, sp, period, nlppform, nlequid)
!
! Setup a spline
!
INTEGER, INTENT(in) :: p, ngauss
DOUBLE PRECISION, INTENT(in) :: grid(:)
TYPE(spline1d), INTENT(out) :: sp
LOGICAL, OPTIONAL, INTENT(in) :: period
LOGICAL, OPTIONAL, INTENT(in) :: nlppform
LOGICAL, OPTIONAL, INTENT(in) :: nlequid
!
DOUBLE COMPLEX :: zc
DOUBLE PRECISION :: leng, xp, factinv, h
INTEGER :: order, nints, i, k
DOUBLE PRECISION :: temp(1:p+1,0:p)
!
! Order of splines
order = p+1
sp%order = order
!
! Dimension of spline space
nints = SIZE(grid)-1
sp%nints = nints
sp%period = .FALSE.
IF( PRESENT(period) ) THEN
sp%period = period
sp%lperiod = grid(nints+1) - grid(1)
END IF
sp%dim = nints+p
!
! Use or not PPFORM
sp%nlppform = .TRUE.
IF( PRESENT(nlppform) ) THEN
sp%nlppform = nlppform
END IF
!
! Determine sequence of knots
IF( ASSOCIATED(sp%knots) ) DEALLOCATE(sp%knots)
ALLOCATE( sp%knots(-p:nints+p) )
sp%knots(0:nints) = grid(:)
!
! Is the grid equidistant ?
IF( PRESENT(nlequid) ) THEN
sp%nlequid = nlequid
ELSE
sp%nlequid = is_equid(grid)
END IF
!
! Coarse to fine mesh mapping for non-equidistant mesh
IF(sp%nlequid) THEN
sp%hinv = 1.0d0/(sp%knots(1)-sp%knots(0))
ELSE
IF(ASSOCIATED(sp%fmap)) DEALLOCATE(sp%fmap)
CALL create_fine(sp%knots(0:nints), h, sp%fmap)
sp%hinv = 1.0d0/h
END IF
!
! Extend knots at both sides of given grid points
IF( sp%period ) THEN
leng = sp%knots(nints) - sp%knots(0)
DO i=-1,-p,-1
sp%knots(i) = sp%knots(nints+i) - leng
END DO
DO i=1,p
sp%knots(nints+i) = sp%knots(i) + leng
END DO
!!$ sp%knots(-p:-1) = sp%knots(nints-p:nints-1) - leng
!!$ sp%knots(nints+1:nints+p) = leng + sp%knots(1:p)
ELSE
sp%knots(-p:-1) = sp%knots(0)
sp%knots(nints+1:nints+p) = sp%knots(nints)
END IF
!
! Precalculated values of all splines and their derivatives at left boundaries
IF( ASSOCIATED(sp%val0) ) DEALLOCATE(sp%val0)
ALLOCATE( sp%val0(0:p, p+1, 1:nints) )
sp%val0 = 0.0d0
DO i=1,nints
xp = sp%knots(i-1) + EPSILON(1.0d0)*ABS(sp%knots(i-1))
CALL basfun_recur(xp, sp, temp, i)
sp%val0(:,:,i) = TRANSPOSE(temp)
END DO
!
factinv = 1.0d0
DO k=2,p ! Divide by k! for use in PPFORM_ALT
factinv = factinv/k
sp%val0(k,:,:) = sp%val0(k,:,:)*factinv
END DO
!
! Case of periodic equidistant splines (translational invariance)
IF(sp%period .AND. sp%nlequid) THEN
IF( ASSOCIATED(sp%valc) ) DEALLOCATE(sp%valc)
ALLOCATE(sp%valc(0:p, p+1))
sp%valc = sp%val0(:,:,1)
END IF
!
! Gauss abscissas and weights
IF( ngauss .GT. 0 ) THEN
IF( ASSOCIATED(sp%gausx) ) DEALLOCATE(sp%gausx)
IF( ASSOCIATED(sp%gausw) ) DEALLOCATE(sp%gausw)
ALLOCATE(sp%gausx(ngauss,nints))
ALLOCATE(sp%gausw(ngauss,nints))
DO i=1,nints
CALL gauleg(sp%knots(i-1), sp%knots(i), &
& sp%gausx(1:ngauss,i), sp%gausw(1:ngauss,i), ngauss)
END DO
END IF
!
! Compute integral of each splines
IF( ASSOCIATED(sp%intspl) ) DEALLOCATE(sp%intspl)
ALLOCATE(sp%intspl(0:sp%dim-1))
CALL calc_integ(sp, sp%intspl)
!
END SUBROUTINE set_spline1d
!===========================================================================
SUBROUTINE init_dft(sp, kmin, kmax, couplings)
!
! Initialize DFT
!
TYPE(spline1d) :: sp
INTEGER, INTENT(in) :: kmin, kmax
INTEGER, INTENT(in), OPTIONAL :: couplings(:)
!
INTEGER :: n, p, dk, k, j, nc
DOUBLE COMPLEX :: zc
!
n = sp%nints
p = sp%order-1
dk = kmax-kmin+1
!
! Check that -N/2 .LE. Kmin .LE. Kmax .LT. N/2
!
IF(kmin.GT.kmax .OR. kmin.LT.-n/2 .OR. kmax.GE.n/2) THEN
WRITE(*,'(a,2i6,a)') 'kmin, kmax =', kmin, kmax, ' erroneous!'
STOP
END IF
!
! The Fourier window
!
sp%dft%n = n
sp%dft%kmin = kmin
sp%dft%kmax = kmax
sp%dft%dk = dk
sp%dft%dx = sp%knots(1) - sp%knots(0)
!
! Precalculate the DFT coefs exp( i(2*pi/N)jk ), k=kmin,kmax, j=0,p
!
IF( ASSOCIATED(sp%dft%coefs) ) DEALLOCATE(sp%dft%coefs)
ALLOCATE(sp%dft%coefs(kmin:kmax,0:p))
zc = EXP( CMPLX(0.0d0, 2.0d0*pi/REAL(n,8),8) )
sp%dft%coefs(:,0) = 1.0d0 ! j=0
sp%dft%coefs(kmin,1) = zc**kmin ! j=1
DO k=kmin+1,kmax
sp%dft%coefs(k,1) = sp%dft%coefs(k-1,1)*zc
END DO
DO j=2,p
sp%dft%coefs(:,j) = sp%dft%coefs(:,1)*sp%dft%coefs(:,j-1)
END DO
!
! Mode couplings: by default use the whole window
!
nc = dk
IF( PRESENT(couplings)) nc = SIZE(couplings)
!
IF(ASSOCIATED(sp%dft%mode_couplings)) DEALLOCATE(sp%dft%mode_couplings)
ALLOCATE(sp%dft%mode_couplings(nc))
!
IF(PRESENT(couplings)) THEN
sp%dft%mode_couplings = couplings
ELSE
sp%dft%mode_couplings = (/ (k,k=kmin,kmax) /)
END IF
END SUBROUTINE init_dft
!===========================================================================
SUBROUTINE set_spline2d(p, ngauss, grid1, grid2, sp, period, nlppform,&
& nlequid)
!
! Setup a 2d spline
!
INTEGER, INTENT(in) :: p(2), ngauss(2)
DOUBLE PRECISION, INTENT(in) :: grid1(:)
DOUBLE PRECISION, INTENT(in) :: grid2(:)
TYPE(spline2d), INTENT(out) :: sp
LOGICAL, OPTIONAL, INTENT(in) :: period(2)
LOGICAL, OPTIONAL, INTENT(in) :: nlppform
LOGICAL, OPTIONAL, INTENT(in) :: nlequid(2)
!
IF(PRESENT(period).AND.PRESENT(nlppform).AND.PRESENT(nlequid)) THEN
CALL set_spline1d(p(1), ngauss(1), grid1, sp%sp1, period(1), nlppform, &
& nlequid(1))
CALL set_spline1d(p(2), ngauss(2), grid2, sp%sp2, period(2), nlppform, &
& nlequid(2))
ELSE IF(PRESENT(period).AND.PRESENT(nlppform)) THEN
CALL set_spline1d(p(1), ngauss(1), grid1, sp%sp1, period(1), nlppform)
CALL set_spline1d(p(2), ngauss(2), grid2, sp%sp2, period(2), nlppform)
ELSE IF(PRESENT(period).AND.PRESENT(nlequid)) THEN
CALL set_spline1d(p(1), ngauss(1), grid1, sp%sp1, period=period(1), nlequid=nlequid(1))
CALL set_spline1d(p(2), ngauss(2), grid2, sp%sp2, period=period(2), nlequid=nlequid(2))
ELSE IF(PRESENT(nlppform).AND.PRESENT(nlequid)) THEN
CALL set_spline1d(p(1), ngauss(1), grid1, sp%sp1, nlppform=nlppform, nlequid=nlequid(1))
CALL set_spline1d(p(2), ngauss(2), grid2, sp%sp2, nlppform=nlppform, nlequid=nlequid(2))
ELSE IF(PRESENT(period)) THEN
CALL set_spline1d(p(1), ngauss(1), grid1, sp%sp1, period(1))
CALL set_spline1d(p(2), ngauss(2), grid2, sp%sp2, period(2))
ELSE IF(PRESENT(nlppform)) THEN
CALL set_spline1d(p(1), ngauss(1), grid1, sp%sp1, nlppform=nlppform)
CALL set_spline1d(p(2), ngauss(2), grid2, sp%sp2, nlppform=nlppform)
ELSE IF(PRESENT(nlequid)) THEN
CALL set_spline1d(p(1), ngauss(1), grid1, sp%sp1, nlequid=nlequid(1))
CALL set_spline1d(p(2), ngauss(2), grid2, sp%sp2, nlequid=nlequid(2))
ELSE
CALL set_spline1d(p(1), ngauss(1), grid1, sp%sp1)
CALL set_spline1d(p(2), ngauss(2), grid2, sp%sp2)
END IF
END SUBROUTINE set_spline2d
!===========================================================================
SUBROUTINE set_spline2d1d(p, ngauss, grid1, grid2, grid3, sp, period, &
& nlppform, nlequid)
!
! Setup a 2d1d spline (for axisymmetric problems)
!
INTEGER, INTENT(in) :: p(3), ngauss(3)
DOUBLE PRECISION, INTENT(in) :: grid1(:)
DOUBLE PRECISION, INTENT(in) :: grid2(:)
DOUBLE PRECISION, INTENT(in) :: grid3(:)
TYPE(spline2d1d), INTENT(out) :: sp
LOGICAL, OPTIONAL, INTENT(in) :: period(3)
LOGICAL, OPTIONAL, INTENT(in) :: nlppform
LOGICAL, OPTIONAL, INTENT(in) :: nlequid(3)
!
IF(PRESENT(period).AND.PRESENT(nlppform).AND.PRESENT(nlequid)) THEN
CALL set_spline2d(p(1:2), ngauss(1:2), grid1, grid2, sp%sp12, period(1:2),&
& nlppform, nlequid(1:2))
CALL set_spline1d(p(3), ngauss(3), grid3, sp%sp3, period(3), nlppform, &
& nlequid(3))
ELSE IF(PRESENT(period).AND.PRESENT(nlppform)) THEN
CALL set_spline2d(p(1:2), ngauss(1:2), grid1, grid2, sp%sp12, period(1:2),&
& nlppform)
CALL set_spline1d(p(3), ngauss(3), grid3, sp%sp3, period(3), nlppform)
ELSE IF(PRESENT(period).AND.PRESENT(nlequid)) THEN
CALL set_spline2d(p(1:2), ngauss(1:2), grid1, grid2, sp%sp12, period=period(1:2),&
& nlequid=nlequid(1:2))
CALL set_spline1d(p(3), ngauss(3), grid3, sp%sp3, period=period(3), &
& nlequid=nlequid(3))
ELSE IF(PRESENT(nlppform).AND.PRESENT(nlequid)) THEN
CALL set_spline2d(p(1:2), ngauss(1:2), grid1, grid2, sp%sp12, nlppform=nlppform,&
& nlequid=nlequid(1:2))
CALL set_spline1d(p(3), ngauss(3), grid3, sp%sp3, nlppform=nlppform, &
& nlequid=nlequid(3))
ELSE IF(PRESENT(period)) THEN
CALL set_spline2d(p(1:2), ngauss(1:2), grid1, grid2, sp%sp12, period(1:2))
CALL set_spline1d(p(3), ngauss(3), grid3, sp%sp3, period(3))
ELSE IF(PRESENT(nlppform)) THEN
CALL set_spline2d(p(1:2), ngauss(1:2), grid1, grid2, sp%sp12, nlppform=nlppform)
CALL set_spline1d(p(3), ngauss(3), grid3, sp%sp3, nlppform=nlppform)
ELSE IF(PRESENT(nlequid)) THEN
CALL set_spline2d(p(1:2), ngauss(1:2), grid1, grid2, sp%sp12, nlequid=nlequid(1:2))
CALL set_spline1d(p(3), ngauss(3), grid3, sp%sp3, nlequid=nlequid(3))
ELSE
CALL set_spline2d(p(1:2), ngauss(1:2), grid1, grid2, sp%sp12)
CALL set_spline1d(p(3), ngauss(3), grid3, sp%sp3)
END IF
END SUBROUTINE set_spline2d1d
!===========================================================================
SUBROUTINE get_dim1(sp, dim, nx, nidbas)
!
! Return spline dimension of 1d spline sp and optionally
! number of knot intervals nx and spline degree nidbas
!
TYPE(spline1d), INTENT(in) :: sp
INTEGER, INTENT(out) :: dim
INTEGER, OPTIONAL, INTENT(out) :: nx, nidbas
dim = sp%dim
IF( PRESENT(nx) ) nx = sp%nints
IF( PRESENT(nidbas) ) nidbas = sp%order - 1
END SUBROUTINE get_dim1
!===========================================================================
SUBROUTINE get_dim2(sp, dim, nx, nidbas)
!
! Return spline dimension of 2d spline sp and optionally
! number of knot intervals nx and spline degree nidbas
!
TYPE(spline2d), INTENT(in) :: sp
INTEGER, INTENT(out) :: dim(2)
INTEGER, OPTIONAL, INTENT(out) :: nx(2), nidbas(2)
!
dim(1) = sp%sp1%dim
IF( PRESENT(nx) ) nx(1) = sp%sp1%nints
IF( PRESENT(nidbas) ) nidbas(1) = sp%sp1%order - 1
!
dim(2) = sp%sp2%dim
IF( PRESENT(nx) ) nx(2) = sp%sp2%nints
IF( PRESENT(nidbas) ) nidbas(2) = sp%sp2%order - 1
END SUBROUTINE get_dim2
!===========================================================================
SUBROUTINE get_gauss(sp, n, i, x, w)
!
! Get Gauss points and weights from spline sp
!
TYPE(spline1d), INTENT(in) :: sp
INTEGER, INTENT(out) :: n
INTEGER, INTENT(in), OPTIONAL :: i
DOUBLE PRECISION, DIMENSION(:), OPTIONAL, INTENT(out) :: x, w
!
n = SIZE(sp%gausx, 1)
IF( PRESENT(i) ) THEN
x(:) = sp%gausx(:,i)
w(:) = sp%gausw(:,i)
END IF
END SUBROUTINE get_gauss
!===========================================================================
SUBROUTINE def_basfun0(xp, sp, fun, left)
!
! Define the basis function and its derivatives at x
! fun(i,j) = (j-1)th derivative of ith basis function.
!
DOUBLE PRECISION, INTENT(in) :: xp
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(out) :: fun(:,:)
INTEGER, OPTIONAL, INTENT(out) :: left
DOUBLE PRECISION :: x
INTEGER :: p, n, kleft
INTEGER :: ierr, j, k
!
CALL locintv(sp, xp, kleft)
CALL basfun(xp, sp, fun, kleft+1)
IF(PRESENT(left)) THEN
left = kleft
END IF
END SUBROUTINE def_basfun0
!===========================================================================
SUBROUTINE basfun0(xp, sp, f, left)
!
! Define the basis function and its derivatives at x, in interval
! [left,left+1], using PPFORM defined by sp%val0., left=1,..,nints
! f(i,j) = jth derivative of ith basis function.
!
DOUBLE PRECISION, INTENT(in) :: xp
+ TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(out) :: f(:,0:)
INTEGER, INTENT(in) :: left ! =1,2,...,nints
- TYPE(spline1d) :: sp
!
INTEGER :: p, n, jdermx, i, jder
DOUBLE PRECISION :: x, h
!
p = sp%order - 1
n = sp%nints
jdermx = SIZE(f,2)-1
!
h = xp-sp%knots(left-1) ! knots are numbered from 0
!
IF(sp%period .AND. sp%nlequid) THEN
DO jder=0,jdermx ! Derivative jder
CALL my_ppval(p, h, sp%valc, jder, f(:,jder))
END DO
ELSE
DO jder=0,jdermx ! Derivative jder
CALL my_ppval(p, h, sp%val0(:,:,left), jder, f(:,jder))
END DO
END IF
CONTAINS
SUBROUTINE my_ppval(p, x, ppform, jder, f)
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x, ppform(:,:)
INTEGER, INTENT(in) :: jder
DOUBLE PRECISION, INTENT(out) :: f(:)
DOUBLE PRECISION :: fact
INTEGER :: j
SELECT CASE (jder)
CASE(0) ! function value
SELECT CASE(p)
CASE(1)
f(1) = ppform(1,1) + x*ppform(2,1)
CASE(2)
f(1) = ppform(1,1) + x*(ppform(2,1)+x*ppform(3,1))
f(2) = ppform(1,2) + x*(ppform(2,2)+x*ppform(3,2))
CASE(3)
f(1) = ppform(1,1) + x*(ppform(2,1)+x*(ppform(3,1)+x*ppform(4,1)))
f(2) = ppform(1,2) + x*(ppform(2,2)+x*(ppform(3,2)+x*ppform(4,2)))
f(3) = ppform(1,3) + x*(ppform(2,3)+x*(ppform(3,3)+x*ppform(4,3)))
CASE(4:)
f(1:p) = ppform(p+1,1:p)
DO j=p,1,-1
f(1:p) = f(1:p)*x + ppform(j,1:p)
END DO
END SELECT
f(p+1) = 1.0d0 - SUM(f(1:p))
CASE(1) ! 1st derivative
SELECT CASE(p)
CASE(1)
f(1) = ppform(2,1)
CASE(2)
f(1) = ppform(2,1) + x*2.d0*ppform(3,1)
f(2) = ppform(2,2) + x*2.d0*ppform(3,2)
CASE(3)
f(1) = ppform(2,1) + x*(2.d0*ppform(3,1)+x*3.0d0*ppform(4,1))
f(2) = ppform(2,2) + x*(2.d0*ppform(3,2)+x*3.0d0*ppform(4,2))
f(3) = ppform(2,3) + x*(2.d0*ppform(3,3)+x*3.0d0*ppform(4,3))
CASE(4:)
f(1:p) = p*ppform(p+1,1:p)
DO j=p-1,1,-1
f(1:p) = f(1:p)*x + j*ppform(j+1,1:p)
END DO
END SELECT
f(p+1) = -SUM(f(1:p))
CASE default ! 2nd and higher derivatives
fact = p-jder
f(1:p) = ppform(p+1,1:p)
DO j=p,jder+1,-1
f(1:p) = f(1:p)/fact*j*x + ppform(j,1:p)
fact = fact-1.0d0
END DO
DO j=2,jder
f(1:p) = f(1:p)*j
END DO
f(p+1) = -SUM(f(1:p))
END SELECT
END SUBROUTINE my_ppval
END SUBROUTINE basfun0
!===========================================================================
SUBROUTINE basfun1(xp, sp, f, left)
!
! Define the basis function and its derivatives at x, in interval i=1,2,
! using PPFORM defined by sp%val0.
! f(i,j,p) = jth derivative of ith basis function at coordinate xp
!
- DOUBLE precision, INTENT(in) :: xp(:)
+ DOUBLE PRECISION, INTENT(in) :: xp(:)
+ TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(out) :: f(0:,0:,:)
INTEGER, INTENT(in) :: left(:) ! =1,2,...,nints
- TYPE(spline1d) :: sp
!
INTEGER :: p, n, kleft, i, j, jder, ierr
INTEGER :: npt, jdermx
DOUBLE PRECISION :: h(SIZE(xp)), temp(SIZE(xp))
DOUBLE PRECISION :: ppform(SIZE(xp),sp%order)
!
p = sp%order - 1
n = sp%nints
npt = SIZE(xp)
jdermx = SIZE(f,2)-1
!
h = xp - sp%knots(left-1) ! knots are numbered from 0
!
IF( sp%period .AND. sp%nlequid) THEN
DO jder=0,jdermx
CALL my_ppval_same(p, h, sp%valc, jder, f(:,jder,1:npt))
END DO
ELSE
DO i=0,p ! Spline i
DO j=1,npt
ppform(j,:) = sp%val0(:,i+1,left(j))
END DO
DO jder=0,jdermx ! Derivative jder
CALL my_ppval(p, h, ppform, jder, temp)
f(i,jder,1:npt) = temp
END DO
END DO
END IF
CONTAINS
!+++
SUBROUTINE my_ppval(p, x, ppform, jder, f)
!
! Compute function and derivatives from the PP representation
! for many points x(:)
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION, INTENT(in) :: ppform(:,:)
INTEGER, INTENT(in) :: jder
DOUBLE PRECISION, INTENT(out) :: f(:)
DOUBLE PRECISION :: fact
INTEGER :: j
SELECT CASE (jder)
CASE(0) ! function value
SELECT CASE(p)
CASE(1)
f(:) = ppform(:,1) + x(:)*ppform(:,2)
CASE(2)
f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*ppform(:,3))
CASE(3)
f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*(ppform(:,3)+x(:)*ppform(:,4)))
CASE(4:)
f(:) = ppform(:,p+1)
DO j=p,1,-1
f(:) = ppform(:,j) + f(:)*x(:)
END DO
END SELECT
CASE(1) ! 1st derivative
SELECT CASE(p)
CASE(1)
f(:) = ppform(:,2)
CASE(2)
f(:) = ppform(:,2) + x(:)*2.d0*ppform(:,3)
CASE(3)
f(:) = ppform(:,2) + x(:)*(2.d0*ppform(:,3)+x(:)*3.0d0*ppform(:,4))
CASE(4:)
f(:) = p*ppform(:,p+1)
DO j=p-1,1,-1
f(:) = f(:)*x(:) + j*ppform(:,j+1)
END DO
END SELECT
CASE default ! 2nd and higher derivatives
f(:) = ppform(:,p+1)
fact = p-jder
DO j=p,jder+1,-1
f(:) = f(:)/fact*j*x(:) + ppform(:,j)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:) = f(:)*j
END DO
END SELECT
END SUBROUTINE my_ppval
!+++
SUBROUTINE my_ppval_same(p, x, ppform, jder, f)
!
! Compute function and derivatives from the PP representation
! for many points x(:), same ppform (translationnal invariant spline)
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION, INTENT(in) :: ppform(:,:)
INTEGER, INTENT(in) :: jder
DOUBLE PRECISION, INTENT(out) :: f(:,:)
DOUBLE PRECISION :: fact
INTEGER :: j,k
SELECT CASE (jder)
!
! function value
CASE(0)
SELECT CASE(p)
CASE(1)
f(1,:) = ppform(1,1) + x(:)*ppform(2,1)
CASE(2)
f(1,:) = ppform(1,1) + x(:)*(ppform(2,1)+x(:)*ppform(3,1))
f(2,:) = ppform(1,2) + x(:)*(ppform(2,2)+x(:)*ppform(3,2))
CASE(3)
f(1,:) = ppform(1,1) + x(:)*(ppform(2,1)+x(:)*(ppform(3,1)+x(:)*ppform(4,1)))
f(2,:) = ppform(1,2) + x(:)*(ppform(2,2)+x(:)*(ppform(3,2)+x(:)*ppform(4,2)))
f(3,:) = ppform(1,3) + x(:)*(ppform(2,3)+x(:)*(ppform(3,3)+x(:)*ppform(4,3)))
CASE(4:)
DO k=1,p
f(k,:) = ppform(p+1,k)
DO j=p,1,-1
f(k,:) = ppform(j,k) + f(k,:)*x(:)
END DO
END DO
END SELECT
f(p+1,:) = 1.0d0 - SUM(f(1:p,:),DIM=1)
!
! 1st derivative
CASE(1)
SELECT CASE(p)
CASE(1)
f(1,:) = ppform(2,1)
CASE(2)
f(1,:) = ppform(2,1) + x(:)*2.d0*ppform(3,1)
f(2,:) = ppform(2,2) + x(:)*2.d0*ppform(3,2)
CASE(3)
f(1,:) = ppform(2,1) + x(:)*(2.d0*ppform(3,1)+x(:)*3.0d0*ppform(4,1))
f(2,:) = ppform(2,2) + x(:)*(2.d0*ppform(3,2)+x(:)*3.0d0*ppform(4,2))
f(3,:) = ppform(2,3) + x(:)*(2.d0*ppform(3,3)+x(:)*3.0d0*ppform(4,3))
CASE(4:)
DO k=1,p
f(k,:) = p*ppform(p+1,k)
DO j=p-1,1,-1
f(k,:) = f(k,:)*x(:) + j*ppform(j+1,k)
END DO
END DO
END SELECT
f(p+1,:) = -SUM(f(1:p,:),DIM=1)
!
! 2nd and higher derivatives
CASE(2:)
DO k=1,p
f(k,:) = ppform(p+1,k)
fact = p-jder
DO j=p,jder+1,-1
f(k,:) = f(k,:)/fact*j*x(:) + ppform(j,k)
fact = fact-1.0d0
END DO
DO j=2,jder
f(k,:) = f(k,:)*j
END DO
END DO
f(p+1,:) = -SUM(f(1:p,:),DIM=1)
END SELECT
END SUBROUTINE my_ppval_same
!+++
END SUBROUTINE basfun1
!===========================================================================
SUBROUTINE def_basfun1(xp, sp, fun, left)
!
! Define the basis function and its derivatives at x
! fun(i,j) = (j-1)th derivative of ith basis function.
!
DOUBLE PRECISION, INTENT(in) :: xp(:)
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(out) :: fun(:,:,:)
INTEGER, OPTIONAL, INTENT(out) :: left(:)
DOUBLE PRECISION :: x(SIZE(xp))
INTEGER :: kleft(SIZE(xp))
INTEGER :: p, n
INTEGER :: ierr, j, k
!
CALL locintv(sp, xp, kleft)
CALL basfun(xp, sp, fun, kleft+1)
IF(PRESENT(left)) THEN
left = kleft
END IF
END SUBROUTINE def_basfun1
!===========================================================================
SUBROUTINE ft_basfun0(xp, sp, ft_f, left)
!
! DFT of basis functions: ft_f(k,j), k=sp%dft%kmin, sp$dft%kmax (modes)
! j=0, p-1 (order of derivative))
!
DOUBLE PRECISION, INTENT(in) :: xp
DOUBLE COMPLEX, INTENT(out) :: ft_f(:,:)
INTEGER, INTENT(in) :: left
TYPE(spline1d) :: sp
DOUBLE PRECISION :: f(sp%order,SIZE(ft_f,2))
!
! Construct all splines on interval [left,left+1] at coordinate xp
CALL basfun(xp, sp, f, left)
!
! DFT of splines
ft_f = MATMUL(sp%dft%coefs, f)
END SUBROUTINE ft_basfun0
!===========================================================================
SUBROUTINE ft_basfun1(xp, sp, ft_f, left)
!
! DFT of basis functions: ft_f(k,j), k=sp%dft%kmin, sp$dft%kmax (modes)
! j=0, p-1 (order of derivative))
! at xp(i)
!
DOUBLE PRECISION, INTENT(in) :: xp(:)
DOUBLE COMPLEX, INTENT(out) :: ft_f(:,:,:)
INTEGER, INTENT(in) :: left(:)
TYPE(spline1d) :: sp
!
INTEGER :: i, n3
DOUBLE PRECISION :: f(sp%order,SIZE(ft_f,2),SIZE(ft_f,3))
!
! Construct all splines on interval [left,left+1] at coordinate xp
CALL basfun(xp, sp, f, left)
!
! DFT of splines
n3 = SIZE(xp)
DO i=1,n3
ft_f(:,:,i) = MATMUL(sp%dft%coefs, f(:,:,i))
END DO
END SUBROUTINE ft_basfun1
!===========================================================================
SUBROUTINE basfun_recur(xp, sp, fun, left)
!
! Define the basis function and its derivatives at x, in interval i=1,2,
! using recurrence construct in function BVALUE
! fun(i,j) = (j-1)th derivative of ith basis function.
!
DOUBLE PRECISION, INTENT(in) :: xp
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(out) :: fun(:,:)
INTEGER, INTENT(in) :: left
DOUBLE PRECISION :: bcoef(1)=1.0d0, bvalue, x
INTEGER :: p, n, kleft
INTEGER :: ierr, j, k
!
p = sp%order - 1
n = sp%nints
fun = 0.0d0
!
IF( sp%period ) THEN ! ** Applly periodicity **
x =sp%knots(0) + MODULO(xp-sp%knots(0), sp%lperiod)
ELSE
x = xp
END IF
!
kleft = left-1
DO j=kleft-p, kleft
DO k=0, SIZE(fun,2)-1
fun(j-kleft+p+1, k+1) = bvalue(sp%knots(j), bcoef, 1, p+1, x, k)
END DO
END DO
END SUBROUTINE basfun_recur
!===========================================================================
SUBROUTINE gauleg(x1,x2,x,w,n)
!
! Compute Gauss-Legendre abscissas and weights in interval [x1, x2]
!
INTEGER, INTENT(in) :: n
DOUBLE PRECISION, INTENT(in) :: x1,x2
DOUBLE PRECISION, INTENT(out) :: x(n),w(n)
DOUBLE PRECISION :: EPS
INTEGER i,j,m
DOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1
!
eps=EPSILON(eps)
m=(n+1)/2
xm=0.5d0*(x2+x1)
xl=0.5d0*(x2-x1)
DO i=1,m
z=COS(3.141592654d0*(i-.25d0)/(n+.5d0))
DO
p1=1.d0; p2=0.d0
DO j=1,n
p3=p2; p2=p1
p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
END DO
pp=n*(z*p1-p2)/(z*z-1.d0)
z1=z
z=z1-p1/pp
IF( ABS(z-z1) .LE. EPS ) EXIT
END DO
x(i)=xm-xl*z
x(n+1-i)=xm+xl*z
w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
w(n+1-i)=w(i)
END DO
END SUBROUTINE gauleg
!===========================================================================
SUBROUTINE gridval1dz(sp, xp, f, jder, c, ppformz)
!
! Compute values or jder-th dervivative of f(x) from ppform
! of spline sp. Recompute the ppform if the optional spline
! coefficients are given.
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: f
INTEGER, INTENT(in) :: jder
DOUBLE COMPLEX, DIMENSION(:), OPTIONAL, INTENT(in) :: c
DOUBLE COMPLEX, DIMENSION(:,:), OPTIONAL :: ppformz
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION :: x(SIZE(xp)), h, fact
INTEGER :: order, nints, i, j, nidbas
INTEGER :: leftx(SIZE(xp))
!
order = sp%order
nints = sp%nints
nidbas = order-1
!
! Compute PPFORM/BCOEFS if spline coefs are passed
!
IF (PRESENT(c)) THEN
IF (sp%nlppform) THEN
IF( PRESENT(ppformz) ) THEN
CALL topp0z(sp, c, ppformz)
ELSE
IF( ASSOCIATED(sp%ppformz) ) DEALLOCATE(sp%ppformz)
ALLOCATE(sp%ppformz(order,nints))
CALL topp0z(sp, c, sp%ppformz)
END IF
ELSE
IF( ASSOCIATED(sp%bcoefsc) ) DEALLOCATE(sp%bcoefsc)
ALLOCATE(sp%bcoefsc(SIZE(c)))
sp%bcoefsc = c
END IF
END IF
!
! Applly periodicity if required
!
IF( sp%period ) THEN
x =sp%knots(0) + MODULO(xp-sp%knots(0), sp%lperiod)
ELSE
x = xp
END IF
!
! Locate the intervals containing x
!
CALL locintv(sp, x, leftx)
!
! Compute function/derivatives
!
IF( sp%nlppform ) THEN ! using PP form
DO i=1,SIZE(x)
IF( PRESENT(ppformz) ) THEN
CALL ppval(sp, x(i), ppformz(:,leftx(i)+1), leftx(i), jder, f(i))
ELSE
CALL ppval(sp, x(i), sp%ppformz(:,leftx(i)+1), leftx(i), jder, f(i))
END IF
END DO
ELSE ! using spline expansion
ALLOCATE(fun(0:nidbas,0:jder))
f = 0.0d0
DO i=1,SIZE(x)
CALL basfun(x(i), sp, fun, leftx(i)+1)
DO j=0,nidbas
f(i) = f(i) + sp%bcoefsc(leftx(i)+j+1)*fun(j,jder)
END DO
END DO
DEALLOCATE(fun)
END IF
!
END SUBROUTINE gridval1dz
!===========================================================================
SUBROUTINE gridval1d(sp, xp, f, jder, c, ppform)
!
! Compute values or jder-th dervivative of f(x) from ppform
! of spline sp. Recompute the ppform if the optional spline
! coefficients are given.
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: f
INTEGER, INTENT(in) :: jder
DOUBLE PRECISION, DIMENSION(:), OPTIONAL, INTENT(in) :: c
DOUBLE PRECISION, DIMENSION(:,:), OPTIONAL :: ppform
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION :: x(SIZE(xp)), h, fact
INTEGER :: order, nints, i, j, nidbas
INTEGER :: leftx(SIZE(xp))
!
order = sp%order
nints = sp%nints
nidbas = order-1
!
! Compute PPFORM/BCOEFS if spline coefs are passed
!
IF (PRESENT(c)) THEN
IF (sp%nlppform) THEN
IF( PRESENT(ppform) ) THEN
CALL topp0(sp, c, ppform)
ELSE
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
ALLOCATE(sp%ppform(order,nints))
CALL topp0(sp, c, sp%ppform)
END IF
ELSE
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
ALLOCATE(sp%bcoefs(SIZE(c)))
sp%bcoefs = c
END IF
END IF
!
! Applly periodicity if required
!
IF( sp%period ) THEN
x =sp%knots(0) + MODULO(xp-sp%knots(0), sp%lperiod)
ELSE
x = xp
END IF
!
! Locate the intervals containing x
!
CALL locintv(sp, x, leftx)
!
! Compute function/derivatives
!
IF( sp%nlppform ) THEN ! using PP form
DO i=1,SIZE(x)
IF( PRESENT(ppform) ) THEN
CALL ppval(sp, x(i), ppform(:,leftx(i)+1), leftx(i), jder, f(i))
ELSE
CALL ppval(sp, x(i), sp%ppform(:,leftx(i)+1), leftx(i), jder, f(i))
END IF
END DO
ELSE ! using spline expansion
ALLOCATE(fun(0:nidbas,0:jder))
f = 0.0d0
DO i=1,SIZE(x)
CALL basfun(x(i), sp, fun, leftx(i)+1)
DO j=0,nidbas
f(i) = f(i) + sp%bcoefs(leftx(i)+j+1)*fun(j,jder)
END DO
END DO
DEALLOCATE(fun)
END IF
!
END SUBROUTINE gridval1d
!===========================================================================
SUBROUTINE def_knots(p, xg, knots, period, nlskip)
!
! Define spline knots for interpolating at sites given by xg
!
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: xg(0:)
DOUBLE PRECISION, POINTER :: knots(:)
LOGICAL, OPTIONAL, INTENT(in) :: period, nlskip
LOGICAL :: kperiod, mlskip
INTEGER :: npt, dim, nx, i, ii
!
kperiod=.FALSE.
mlskip = .TRUE.
IF( PRESENT(period) ) kperiod=period
IF( PRESENT(nlskip) ) mlskip = nlskip
!
! Periodic splines
!
IF( kperiod ) THEN
nx = SIZE(xg) -1
IF( ASSOCIATED(knots) ) DEALLOCATE(knots)
ALLOCATE(knots(0:nx))
IF( MODULO(p,2) .NE. 0 ) THEN ! Odd degree
knots(0:nx) = xg(0:nx)
ELSE ! Even degree
DO i=1,nx
knots(i) = 0.5d0*(xg(i-1)+xg(i))
END DO
knots(0) = knots(nx) - (xg(nx)-xg(0))
END IF
RETURN
END IF
!
! Non-periodic splines
!
npt = SIZE(xg)
dim = npt
IF( .NOT. mlskip ) THEN
dim = dim + 2*(p/2) ! Add BC on derivatives
END IF
nx = dim-p
IF( ASSOCIATED(knots) ) DEALLOCATE(knots)
ALLOCATE(knots(0:nx))
!
knots(0) = xg(0)
knots(nx) = xg(npt-1)
!
IF( MODULO(p,2) .EQ. 0 ) THEN
ii = 0
IF( mlskip ) ii = p/2 ! skip first p/2 intervals
DO i=1,nx-1
ii = ii+1
knots(i) = (xg(ii)+xg(ii-1))/2
END DO
ELSE
ii = 0
IF( mlskip ) ii = (p-1)/2 ! skip (p-1)/2 points after the first point ii=0
DO i=1,nx-1
ii = ii+1
knots(i) = xg(ii)
END DO
END IF
!
END SUBROUTINE def_knots
!===========================================================================
SUBROUTINE allsplines(sp, xpt, splines)
!
! Return all splines defined on points xpt
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: xpt(:)
DOUBLE PRECISION, POINTER :: splines(:,:)
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
INTEGER :: i, n, left, dim, p
!
p = sp%order - 1
dim = sp%dim
n = SIZE(xpt)
!
IF( ASSOCIATED(splines) ) DEALLOCATE(splines)
ALLOCATE(splines(n,dim), fun(p+1,1))
splines = 0.0d0
DO i=1,n
CALL locintv(sp, xpt(i), left)
CALL basfun(xpt(i), sp, fun, left+1)
splines(i,left+1:left+p+1) = fun(1:p+1,1)
END DO
DEALLOCATE(fun)
END SUBROUTINE allsplines
!===========================================================================
SUBROUTINE set_splcoef1d(p, x, sp, period, ibc)
!
! Setup 1d interpolation matrix for spline of degree p
!
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x(:)
TYPE(spline1d), INTENT(out) :: sp
LOGICAL, OPTIONAL, INTENT(in) :: period
INTEGER, OPTIONAL :: ibc(:,:)
!
LOGICAL :: kperiod
!
kperiod = .FALSE.
IF( PRESENT(period) ) kperiod = period
!
IF( kperiod ) THEN
CALL splcoefp_setup(p, x, sp)
ELSE
IF( PRESENT(ibc) ) THEN
CALL splcoef_setup(p, x, sp, ibc)
ELSE
CALL splcoef_setup(p, x, sp)
END IF
END IF
END SUBROUTINE set_splcoef1d
!===========================================================================
SUBROUTINE set_splcoef2d(p, x1, x2, sp, period, ibc1, ibc2)
!
! Setup 2d interpolation matrix for spline of degree p
!
INTEGER, INTENT(in) :: p(2)
DOUBLE PRECISION, INTENT(in) :: x1(:), x2(:)
TYPE(spline2d), INTENT(out) :: sp
LOGICAL, OPTIONAL, INTENT(in) :: period(2)
INTEGER, OPTIONAL :: ibc1(:,:),ibc2(:,:)
!
LOGICAL :: kperiod(2)
!
kperiod = .FALSE.
IF( PRESENT(period) ) kperiod = period
!
! Direction 1
IF( kperiod(1) ) THEN
CALL splcoefp_setup(p(1), x1, sp%sp1)
ELSE
IF( PRESENT(ibc1) ) THEN
CALL splcoef_setup(p(1), x1, sp%sp1, ibc1)
ELSE
CALL splcoef_setup(p(1), x1, sp%sp1)
END IF
END IF
!
! Direction 2
IF( kperiod(2) ) THEN
CALL splcoefp_setup(p(2), x2, sp%sp2)
ELSE
IF( PRESENT(ibc2) ) THEN
CALL splcoef_setup(p(2), x2, sp%sp2, ibc2)
ELSE
CALL splcoef_setup(p(2), x2, sp%sp2)
END IF
END IF
END SUBROUTINE set_splcoef2d
!===========================================================================
SUBROUTINE get_splcoef1(sp, f, c, fbc)
!
! Compute the spline coefficients c from grid values f
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: f(:)
DOUBLE PRECISION, INTENT(out) :: c(:)
DOUBLE PRECISION, INTENT(in), OPTIONAL :: fbc(:,:)
!
IF( sp%period ) THEN
CALL splcoefp1(sp, f, c)
ELSE
IF( PRESENT(fbc) ) THEN
CALL splcoef1(sp, f, c, fbc)
ELSE
CALL splcoef1(sp, f, c)
END IF
END IF
END SUBROUTINE get_splcoef1
!===========================================================================
SUBROUTINE get_splcoef2d(sp, f, c, fbc1, fbc2)
!
! Compute the spline coefficients c from 2d grid values f
!
TYPE(spline2d) :: sp
DOUBLE PRECISION, INTENT(in) :: f(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:,:)
DOUBLE PRECISION :: ctr(SIZE(c,2), SIZE(f,1))
DOUBLE PRECISION, INTENT(in), OPTIONAL :: fbc1(:,:,:), fbc2(:,:,:)
DOUBLE PRECISION, DIMENSION(:, :) , ALLOCATABLE :: c_fbc1_left, c_fbc1_right
DOUBLE PRECISION, DIMENSION(:, :, :), ALLOCATABLE :: c_fbc1_all
!
! Along direction 2
!
IF( PRESENT(fbc2) ) THEN
CALL get_splcoefn(sp%sp2, TRANSPOSE(f), ctr, fbc2)
ELSE
CALL get_splcoefn(sp%sp2, TRANSPOSE(f), ctr)
END IF
!
! Along direction 1
!
IF( PRESENT(fbc1) ) THEN
ALLOCATE( c_fbc1_left(SIZE(c, 2), SIZE(fbc1, 2)))
ALLOCATE(c_fbc1_right(SIZE(c, 2), SIZE(fbc1, 2)))
ALLOCATE(c_fbc1_all(2, SIZE(fbc1, 2), SIZE(c, 2)))
CALL get_splcoefn(sp%sp2, TRANSPOSE(fbc1(1, :, :)), c_fbc1_left )
CALL get_splcoefn(sp%sp2, TRANSPOSE(fbc1(2, :, :)), c_fbc1_right)
c_fbc1_all(1, :, :) = TRANSPOSE(c_fbc1_left )
c_fbc1_all(2, :, :) = TRANSPOSE(c_fbc1_right)
CALL get_splcoefn(sp%sp1, TRANSPOSE(ctr), c, c_fbc1_all)
DEALLOCATE(c_fbc1_left, c_fbc1_right, c_fbc1_all)
ELSE
CALL get_splcoefn(sp%sp1, TRANSPOSE(ctr), c)
END IF
!
END SUBROUTINE get_splcoef2d
!===========================================================================
SUBROUTINE get_splcoef2dz(sp, f, c)
!
! Compute the spline coefficients c from 2d grid values f
!
TYPE(spline2d) :: sp
DOUBLE COMPLEX, INTENT(in) :: f(:,:)
DOUBLE COMPLEX, INTENT(out) :: c(:,:)
DOUBLE PRECISION, DIMENSION(SIZE(c,1), SIZE(c,2),2) :: pc
!
CALL get_splcoef2d(sp, REAL(f), pc(:,:,1))
CALL get_splcoef2d(sp, AIMAG(f), pc(:,:,2))
c(:,:) = CMPLX(pc(:,:,1),pc(:,:,2))
END SUBROUTINE get_splcoef2dz
!===========================================================================
SUBROUTINE get_splcoef1z(sp, f, c, fbc)
!
! Compute the spline coefficients c from grid values f
!
TYPE(spline1d) :: sp
DOUBLE COMPLEX, INTENT(in) :: f(:)
DOUBLE COMPLEX, INTENT(out) :: c(:)
DOUBLE COMPLEX, INTENT(in), OPTIONAL :: fbc(:,:)
DOUBLE PRECISION :: pf(SIZE(f),2), pc(SIZE(c),2)
DOUBLE PRECISION, ALLOCATABLE :: pfbc(:,:,:)
!
pf(:,1) = REAL(f(:))
pf(:,2) = AIMAG(f(:))
IF(PRESENT(fbc)) THEN
ALLOCATE(pfbc(SIZE(fbc,1),SIZE(fbc,2),2))
pfbc(:,:,1) = REAL(fbc(:,:))
pfbc(:,:,2) = AIMAG(fbc(:,:))
CALL get_splcoefn(sp, pf, pc, pfbc)
DEALLOCATE(pfbc)
ELSE
CALL get_splcoefn(sp, pf, pc)
END IF
c(:) = CMPLX(pc(:,1), pc(:,2))
!
END SUBROUTINE get_splcoef1z
!===========================================================================
SUBROUTINE get_splcoefn(sp, f, c, fbc)
!
! Compute the spline coefficients c from grid values f
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: f(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:,:)
DOUBLE PRECISION, INTENT(in), OPTIONAL :: fbc(:,:,:)
!
IF( sp%period ) THEN
CALL splcoefpn(sp, f, c)
ELSE
IF( PRESENT(fbc) ) THEN
CALL splcoefn(sp, f, c, fbc)
ELSE
CALL splcoefn(sp, f, c)
END IF
END IF
END SUBROUTINE get_splcoefn
!===========================================================================
SUBROUTINE splcoef_setup(p, x, sp, ibc)
!
! Setup the interpolation matrix
! for spline of degree p
!
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x(:)
TYPE(spline1d), INTENT(out) :: sp
INTEGER, OPTIONAL :: ibc(:,:)
DOUBLE PRECISION, POINTER :: knots(:)=>NULL(), arow(:)=>NULL(), &
& fun(:,:)=>NULL()
INTEGER :: nx, dim, kl, ku, rank
INTEGER :: i, left, ishift
LOGICAL :: nlskip
!
! Type of Boundary Conditions
nlskip = .TRUE.
ishift = 0
IF( PRESENT(ibc) ) THEN
nlskip = .FALSE.
ishift = p/2
END IF
!
! Set up spline
nx = SIZE(x) - 1 ! X is the interpolation sites
CALL def_knots(p, x, knots, nlskip=nlskip)
CALL set_spline(p, 0, knots, sp)
sp%nsites = nx + 1 ! Store away the number of interpolation sites
DEALLOCATE(knots)
!
! Set up interpolation matrix
dim = sp%dim
kl = MAX(p-1,0)
ku = MAX(p-1,0)
rank = dim
CALL init(kl, ku, rank, 0, sp%mat)
!!$ WRITE(*,'(a,3i6)') 'Interpolation matrix:, kl, ku, rank ', kl, ku, rank
!
! COMPUTE matrix row by row
ALLOCATE(arow(dim), fun(p+1,0:p))
DO i=1,SIZE(x)
arow = 0.0d0
CALL locintv(sp, x(i), left)
CALL basfun(x(i), sp, fun(:,0:0), left+1)
arow(left+1:left+p+1) = fun(1:p+1,0)
CALL putrow(sp%mat, i+ishift, arow)
!!$ WRITE(*,'(i5,13f8.3)') i+ishift, arow
END DO
!
! Add BC if specified
IF( PRESENT(ibc) ) THEN
CALL locintv(sp, x(1), left)
CALL basfun(x(1), sp, fun, left+1) ! BC at the left side
DO i=1,p/2
arow = 0.0d0
arow(left+1:left+p+1) = fun(1:p+1,ibc(1,i))
CALL putrow(sp%mat, i, arow)
!!$ WRITE(*,'(i5,13f8.3)') i, arow
END DO
CALL locintv(sp, x(SIZE(x)), left)
CALL basfun(x(SIZE(x)), sp, fun, left+1) ! BC at the right side
DO i=1,p/2
arow = 0.0d0
arow(left+1:left+p+1) = fun(1:p+1,ibc(2,i))
CALL putrow(sp%mat, dim-i+1, arow)
!!$ WRITE(*,'(i5,13f8.3)') dim-i+1, arow
END DO
END IF
DEALLOCATE(arow, fun)
!
! Factor the matrix
CALL factor(sp%mat)
!
END SUBROUTINE splcoef_setup
!===========================================================================
SUBROUTINE splcoefp_setup(p, x, sp)
!
! Set up the interpolation matrix
! for periodic case
!
USE matrix
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x(:)
TYPE(spline1d), INTENT(out) :: sp
!
TYPE(gemat) :: hmat
DOUBLE PRECISION, POINTER :: knots(:)=>NULL(), arow(:)=>NULL(), &
& fun(:,:)=>NULL()
!!$ DOUBLE PRECISION, POINTER :: arr2d(:,:)=>null()
INTEGER :: nx, kl, ku, rank, mr, nc
INTEGER :: i, left, j, jj
DOUBLE PRECISION, PARAMETER :: zero=0.0d0, one=1.0d0
!________________________________________________________________________________
!
! Set up spline
!
nx = SIZE(x) - 1 ! X is the interpolation sites
CALL def_knots(p, x, knots, period=.TRUE.)
CALL set_spline(p, 0, knots, sp, period=.TRUE.)
sp%nsites = nx + 1 ! Store away the number of interpolation sites
DEALLOCATE(knots)
!________________________________________________________________________________
!
! Set up interpolation matrix sp%matp
!
kl = MAX(p/2,0)
ku = kl
rank = nx
!!$ WRITE(*,'(a,3i6)') 'Interpolation matrix:, kl, ku, rank ', kl, ku, rank
!
CALL init(kl, ku, rank, 0, sp%matp%mat) ! matp%mat is a GB matrix
ALLOCATE(sp%matp%matu(rank, kl+ku), sp%matp%matvt(kl+ku,rank))
!
sp%matp%matu = zero
sp%matp%matvt = zero ! kl = ku = 2
DO j=1,kl ! [ 1 0 0 . . . . ]
sp%matp%matu(j,j) = one ! [ 0 1 0 . . . . ]
END DO ! [ . 0 . . . . . ]
DO j=1,ku ! [ . . . . . 0 . ]
i=rank-ku+j ! [ . . . . 0 1 0 ]
sp%matp%matu(i,kl+j) = one ! [ . . . . 0 0 1 ]
END DO
!________________________________________________________________________________
!
! COMPUTE matrix row by row
!
ALLOCATE(arow(rank), fun(p+1,1))
DO i=1,rank
arow = zero
CALL locintv(sp, x(i), left)
CALL basfun(x(i), sp, fun, left+1)
left = left-p/2
DO j=0,p
jj = MODULO(left+j, rank) + 1
arow(jj) = fun(j+1,1)
END DO
CALL putrow(sp%matp%mat, i, arow)
IF( i .LE. kl ) THEN
sp%matp%matvt(i,rank-kl+1:rank) = arow(rank-kl+1:rank)
ELSE IF ( i .GE. rank-ku+1 ) THEN
j = i-(rank-ku+1) + 1
sp%matp%matvt(kl+j,1:ku) = arow(1:ku)
END IF
!!$ WRITE(8, '(i5, 12(1pe12.3))') left, x(i), fun
!!$ WRITE(*,'(i5, 12(1pe12.3))') i, arow
END DO
DEALLOCATE(arow, fun)
!
!!$ PRINT*, 'Matrix U, V'
!!$ DO i=1,rank
!!$ WRITE(*, '(i5,12(1pe12.3))') i, sp%matp%matu(i,:), sp%matp%matvt(:,i)
!!$ END DO
!!$ ALLOCATE(arr2d(rank,rank))
!!$ arr2d = MATMUL(sp%matp%matu, sp%matp%matvt)
!!$ PRINT*, 'Product U*V^T'
!!$ DO i=1,rank
!!$ WRITE(*, '(i5,12(1pe12.3))') i, arr2d(i,:)
!!$ END DO
!!$ DEALLOCATE(arr2d)
!________________________________________________________________________________
!
! Factorisation
!
! Factor A
CALL factor(sp%matp%mat)
!
! For constant and linear splines, A is diagnonal!
! Should skip the rest
!
IF( kl.EQ.0 .OR. ku.EQ.0 ) THEN
RETURN
END IF
!
! U <-- A^(-1) * U
CALL bsolve(sp%matp%mat, sp%matp%matu)
!
! H <-- 1 + V^T * U
mr = SIZE(sp%matp%matvt, 1)
nc = SIZE(sp%matp%matvt, 2)
CALL init(mr, 0, hmat) ! hmat is initialized to 0!
DO i=1,mr
hmat%val(i,i) = one
END DO
CALL dgemm('N', 'N', mr, mr, nc, one, sp%matp%matvt, mr, &
& sp%matp%matu, nc, one, hmat%val, mr)
!
!!$ hmat%val = MATMUL(sp%matp%matvt, sp%matp%matu)
!!$ DO i=1,kl+ku
!!$ hmat%val(i,i) = 1.0d0 + hmat%val(i,i)
!!$ END DO
!
! V^T <-- H^(-1) V^T
CALL factor(hmat)
CALL bsolve(hmat, sp%matp%matvt)
CALL destroy(hmat)
!
END SUBROUTINE splcoefp_setup
!===========================================================================
SUBROUTINE splcoef1(sp, f, c, fbc)
!
! Compute the spline coefficients c from grid values f and BC fbc
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: f(:)
DOUBLE PRECISION, INTENT(out) :: c(:)
DOUBLE PRECISION, INTENT(in), OPTIONAL :: fbc(:,:)
INTEGER :: p, dim, i, ishift
!
p = sp%order-1
dim = sp%dim
!
! BC at left and right boundary
ishift = 0
IF( PRESENT(fbc) ) THEN
DO i=1,p/2
c(i) = fbc(1,i) ! Left boundary
c(dim-i+1) = fbc(2,i) ! Right boundary
END DO
ishift = p/2
END IF
!
! Interior points
DO i=1,sp%nsites
c(i+ishift) = f(i)
END DO
! WRITE(*,'(a/(13f8.3))') 'RHS', c
!
! Solve for the interpolation coefs. using the factored sp%mat
CALL bsolve(sp%mat, c)
!
END SUBROUTINE splcoef1
!===========================================================================
SUBROUTINE splcoefn(sp, f, c, fbc)
!
! Compute the spline coefficients c from grid values f and BC fbc
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: f(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:,:)
DOUBLE PRECISION, INTENT(in), OPTIONAL :: fbc(:,:,:)
INTEGER :: p, dim, i, ishift
!
p = sp%order-1
dim = sp%dim
!
! BC at left and right boundary
ishift = 0
IF( PRESENT(fbc) ) THEN
DO i=1,p/2
c(i,:) = fbc(1,i,:) ! Left boundary
c(dim-i+1,:) = fbc(2,i,:) ! Right boundary
END DO
ishift = p/2
END IF
!
! Interior points ! c(:,j) for j>SIZE(f,2) could be anything
DO i=1,sp%nsites ! (periodicity in the 2nd dimension)!
c(i+ishift,1:SIZE(f,2)) = f(i,:)
END DO
!!$ WRITE(*,'(a/(13f8.3))') 'RHS', c
!
! Solve for the interpolation coefs. using the factored sp%mat
CALL bsolve(sp%mat, c)
!
END SUBROUTINE splcoefn
!===========================================================================
SUBROUTINE splcoefp1(sp, f, c)
!
! Compute the spline coefficient c from grid values f
! f(x) is periodic
!
USE matrix
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: f(:)
DOUBLE PRECISION, INTENT(out) :: c(:)
!
DOUBLE PRECISION, POINTER :: arow(:), brow(:)
DOUBLE PRECISION, PARAMETER :: zero=0.0d0, one=1.0d0, minus1=-1.0d0
INTEGER :: dim, p, rank, bandw
INTEGER :: i, j
!________________________________________________________________________________
!
p = sp%order-1
rank = sp%nints
bandw = SIZE(sp%matp%matvt,1)
!
! Solve the interpolation system
!
! Solve Ay = f
ALLOCATE(arow(rank), brow(bandw))
!
arow(1:rank) = f(1:rank)
CALL bsolve(sp%matp%mat, arow)
!
! For constant and linear splines, A is diagnonal!
! Should skip the rest
!
IF( p.LE.1 ) GOTO 100
!
!
! t = V^T*y
CALL dgemv('N', bandw, rank, one, sp%matp%matvt, bandw, arow, 1, zero, &
& brow, 1)
!
! y = y - Ut
CALL dgemv('N', rank, bandw, minus1, sp%matp%matu, rank, brow, 1, one, &
& arow, 1)
!
100 CONTINUE
!
! Interpolation coefficients
dim = sp%dim
DO i=1,dim
j = MODULO(i-1-p/2, rank) + 1
c(i) = arow(j)
END DO
!
DEALLOCATE(arow,brow)
!
END SUBROUTINE splcoefp1
!===========================================================================
SUBROUTINE splcoefpn(sp, f, c)
!
! Compute the spline coefficient c from grid values f
! f(x) is periodic
!
USE matrix
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: f(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:,:)
!
DOUBLE PRECISION, POINTER :: arow(:,:), brow(:,:)
DOUBLE PRECISION, PARAMETER :: zero=0.0d0, one=1.0d0, minus1=-1.0d0
INTEGER :: p, dim, rank, bandw, nrhs
INTEGER :: i, j, k
!________________________________________________________________________________
!
p = sp%order-1
rank = sp%nints
bandw = SIZE(sp%matp%matvt,1)
nrhs = SIZE(f,2)
!
!
! Solve the interpolation system
!
! Solve Ay = f
ALLOCATE(arow(rank,nrhs), brow(bandw,nrhs))
!
arow(1:rank,1:nrhs) = f(1:rank,1:nrhs)
CALL bsolve(sp%matp%mat, arow)
!
! For constant and linear splines, A is diagnonal!
! Should skip the rest
!
IF( p.LE.1 ) GOTO 100
!
!
! t = V^T*y
CALL dgemm('N', 'N', bandw, nrhs, rank, one, sp%matp%matvt, bandw, arow, &
& rank, zero, brow, bandw)
!
! y = y - Ut
CALL dgemm('N', 'N', rank, nrhs, bandw, minus1, sp%matp%matu, rank, brow, &
& bandw, one, arow, rank)
!
100 CONTINUE
!
! Interpolation coefficients
dim = sp%dim
DO k=1,nrhs
DO i=1,dim
j = MODULO(i-1-p/2, rank) + 1
c(i,k) = arow(j,k)
END DO
END DO
!
DEALLOCATE(arow,brow)
!
END SUBROUTINE splcoefpn
!
!===========================================================================
SUBROUTINE topp0(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: c(:)
DOUBLE PRECISION, INTENT(out) :: ppform(0:,:)
INTEGER :: p, nints, i, j, k
!
p = sp%order - 1
nints = sp%nints
!
ppform = 0.0d0
DO i=1,nints ! on each knot interval
DO j=1,p+1 ! all spline in interval i
DO k=0,p ! k_th derivatives
ppform(k,i) = ppform(k,i) + sp%val0(k,j,i)*c(j+i-1)
END DO
END DO
END DO
!
END SUBROUTINE topp0
!===========================================================================
SUBROUTINE topp0z(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE COMPLEX, INTENT(in) :: c(:)
DOUBLE COMPLEX, INTENT(out) :: ppform(0:,:)
INTEGER :: p, nints, i, j, k
!
p = sp%order - 1
nints = sp%nints
!
ppform = (0.0d0, 0.0d0)
DO i=1,nints ! on each knot interval
DO j=1,p+1 ! all spline in interval i
DO k=0,p ! k_th derivatives
ppform(k,i) = ppform(k,i) + sp%val0(k,j,i)*c(j+i-1)
END DO
END DO
END DO
!
END SUBROUTINE topp0z
!===========================================================================
SUBROUTINE topp1(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d,:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: c(:,:)
DOUBLE PRECISION, INTENT(out) :: ppform(:,:,:)
INTEGER :: m
!
DO m=1,SIZE(c,2)
CALL topp0(sp, c(:,m), ppform(m,:,:))
END DO
!
END SUBROUTINE topp1
!===========================================================================
SUBROUTINE topp1z(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d,:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE COMPLEX, INTENT(in) :: c(:,:)
DOUBLE COMPLEX, INTENT(out) :: ppform(:,:,:)
INTEGER :: m
!
DO m=1,SIZE(c,2)
CALL topp0z(sp, c(:,m), ppform(m,:,:))
END DO
!
END SUBROUTINE topp1z
!===========================================================================
SUBROUTINE topp2(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d,:,:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: c(:,:,:)
DOUBLE PRECISION, INTENT(out) :: ppform(:,:,:,:)
INTEGER :: m, mm
!
DO mm=1,SIZE(c,3)
DO m=1,SIZE(c,2)
CALL topp0(sp, c(:,m,mm), ppform(m,mm,:,:))
END DO
END DO
!
END SUBROUTINE topp2
!===========================================================================
SUBROUTINE topp2z(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d,:,:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE COMPLEX, INTENT(in) :: c(:,:,:)
DOUBLE COMPLEX, INTENT(out) :: ppform(:,:,:,:)
INTEGER :: m, mm
!
DO mm=1,SIZE(c,3)
DO m=1,SIZE(c,2)
CALL topp0z(sp, c(:,m,mm), ppform(m,mm,:,:))
END DO
END DO
!
END SUBROUTINE topp2z
!===========================================================================
SUBROUTINE topp3(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d,:,:,:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: c(:,:,:,:)
DOUBLE PRECISION, INTENT(out) :: ppform(:,:,:,:,:)
INTEGER :: m, mm, mmm
!
DO mmm=1,SIZE(c,4)
DO mm=1,SIZE(c,3)
DO m=1,SIZE(c,2)
CALL topp0(sp, c(:,m,mm,mmm), ppform(m,mm,mmm,:,:))
END DO
END DO
END DO
!
END SUBROUTINE topp3
!===========================================================================
SUBROUTINE topp3z(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d,:,:,:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE COMPLEX, INTENT(in) :: c(:,:,:,:)
DOUBLE COMPLEX, INTENT(out) :: ppform(:,:,:,:,:)
INTEGER :: m, mm, mmm
!
DO mmm=1,SIZE(c,4)
DO mm=1,SIZE(c,3)
DO m=1,SIZE(c,2)
CALL topp0z(sp, c(:,m,mm,mmm), ppform(m,mm,mmm,:,:))
END DO
END DO
END DO
!
END SUBROUTINE topp3z
!===========================================================================
SUBROUTINE topp4(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d,:,:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: c(:,:,:,:,:)
DOUBLE PRECISION, INTENT(out) :: ppform(:,:,:,:,:,:)
INTEGER :: m, mm, mmm, mmmm
!
DO mmmm=1,SIZE(c,5)
DO mmm=1,SIZE(c,4)
DO mm=1,SIZE(c,3)
DO m=1,SIZE(c,2)
CALL topp0(sp, c(:,m,mm,mmm,mmmm), ppform(m,mm,mmm,mmmm,:,:))
END DO
END DO
END DO
END DO
!
END SUBROUTINE topp4
!===========================================================================
SUBROUTINE topp4z(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d,:,:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE COMPLEX, INTENT(in) :: c(:,:,:,:,:)
DOUBLE COMPLEX, INTENT(out) :: ppform(:,:,:,:,:,:)
INTEGER :: m, mm, mmm, mmmm
!
DO mmmm=1,SIZE(c,5)
DO mmm=1,SIZE(c,4)
DO mm=1,SIZE(c,3)
DO m=1,SIZE(c,2)
CALL topp0z(sp, c(:,m,mm,mmm,mmmm), ppform(m,mm,mmm,mmmm,:,:))
END DO
END DO
END DO
END DO
!
END SUBROUTINE topp4z
!===========================================================================
SUBROUTINE ppval0(sp, x, ppform, left, jder, f)
!
! Compute function and derivatives from the PP representation
! f(x) = \sum_{k=0}^p ppform(k)*(x-t_{left})^k
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x, ppform(:)
INTEGER, INTENT(in) :: left, jder
DOUBLE PRECISION, INTENT(out) :: f
DOUBLE PRECISION :: h, fact
INTEGER :: j, order
!
order = sp%order ! Polynomial degree p + 1
!
h = x-sp%knots(left)
f = 0.0d0
IF( jder .LT. 0 .OR. jder .GE. order ) RETURN
!
SELECT CASE (jder)
CASE(0) ! function value
DO j=order,1,-1
f = f*h + ppform(j)
END DO
CASE(1) ! 1st derivative
DO j=order-1,1,-1
f = f*h + j*ppform(j+1)
END DO
CASE default ! 2nd and higher derivatives
fact = order-jder
DO j=order,jder+1,-1
f = f/fact*j*h + ppform(j)
fact = fact-1.0d0
END DO
DO j=2,jder
f = f*j
END DO
END SELECT
END SUBROUTINE ppval0
!===========================================================================
SUBROUTINE ppval0z(sp, x, ppform, left, jder, f)
!
! Compute function and derivatives from the PP representation
! f(x) = \sum_{k=0}^p ppform(k)*(x-t_{left})^k
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE COMPLEX, INTENT(in) :: ppform(:)
INTEGER, INTENT(in) :: left, jder
DOUBLE COMPLEX, INTENT(out) :: f
DOUBLE PRECISION :: h, fact
INTEGER :: j, order
!
order = sp%order ! Polynomial degree p + 1
!
h = x-sp%knots(left)
f = (0.0d0,0.0d0)
IF( jder .LT. 0 .OR. jder .GE. order ) RETURN
!
SELECT CASE (jder)
CASE(0) ! function value
DO j=order,1,-1
f = f*h + ppform(j)
END DO
CASE(1) ! 1st derivative
DO j=order-1,1,-1
f = f*h + j*ppform(j+1)
END DO
CASE default ! 2nd and higher derivatives
fact = order-jder
DO j=order,jder+1,-1
f = f/fact*j*h + ppform(j)
fact = fact-1.0d0
END DO
DO j=2,jder
f = f*j
END DO
END SELECT
END SUBROUTINE ppval0z
!===========================================================================
SUBROUTINE ppval0z_n(sp, x, ppform, left, jder, f)
!
! PPVAL0Z for many points x(:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE COMPLEX, INTENT(in) :: ppform(:,:)
INTEGER, INTENT(in) :: left(:), jder
DOUBLE COMPLEX, INTENT(out) :: f(:)
DOUBLE PRECISION :: h(SIZE(x)), fact
INTEGER :: j, order
!
order = sp%order ! Polynomial degree p + 1
!
h(:) = x(:)-sp%knots(left(:))
f(:) = 0.0d0
IF( jder .LT. 0 .OR. jder .GE. order ) RETURN
!
SELECT CASE (jder)
CASE(0) ! function value
DO j=order,1,-1
f(:) = f(:)*h(:) + ppform(:,j)
END DO
CASE(1) ! 1st derivative
DO j=order-1,1,-1
f(:) = f(:)*h(:) + j*ppform(:,j+1)
END DO
CASE default ! 2nd and higher derivatives
fact = order-jder
DO j=order,jder+1,-1
f(:) = f(:)/fact*j*h(:) + ppform(:,j)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:) = f(:)*j
END DO
END SELECT
END SUBROUTINE ppval0z_n
!===========================================================================
SUBROUTINE ppval0_n(sp, x, ppform, left, jder, f)
!
! PPVAL0 for many points x(:)
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x(:), ppform(:,:)
INTEGER, INTENT(in) :: left(:), jder
DOUBLE PRECISION, INTENT(out) :: f(:)
DOUBLE PRECISION :: h(SIZE(x)), fact
INTEGER :: j, order
!
order = sp%order ! Polynomial degree p + 1
!
h(:) = x(:)-sp%knots(left(:))
f(:) = 0.0d0
IF( jder .LT. 0 .OR. jder .GE. order ) RETURN
!
SELECT CASE (jder)
CASE(0) ! function value
DO j=order,1,-1
f(:) = f(:)*h(:) + ppform(:,j)
END DO
CASE(1) ! 1st derivative
DO j=order-1,1,-1
f(:) = f(:)*h(:) + j*ppform(:,j+1)
END DO
CASE default ! 2nd and higher derivatives
fact = order-jder
DO j=order,jder+1,-1
f(:) = f(:)/fact*j*h(:) + ppform(:,j)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:) = f(:)*j
END DO
END SELECT
END SUBROUTINE ppval0_n
!===========================================================================
SUBROUTINE ppval1(sp, x, ppform, left, jder, f)
!
! Compute function and derivatives from the PP representation
! f(x) = \sum_{k=0}^p ppform(k)*(x-t_{left})^k
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x, ppform(:,:)
INTEGER, INTENT(in) :: left, jder
DOUBLE PRECISION, INTENT(out) :: f(:)
DOUBLE PRECISION :: h, fact
INTEGER :: j, order
!
order = sp%order ! Polynomial degree p + 1
!
h = x-sp%knots(left)
f = 0.0d0
IF( jder .LT. 0 .OR. jder .GE. order ) RETURN
!
SELECT CASE (jder)
CASE(0) ! function value
DO j=order,1,-1
f(:) = f(:)*h + ppform(j,:)
END DO
CASE(1) ! 1st derivative
DO j=order-1,1,-1
f(:) = f(:)*h + j*ppform(j+1,:)
END DO
CASE default ! 2nd and higher derivatives
fact = order-jder
DO j=order,jder+1,-1
f(:) = f(:)/fact*j*h + ppform(j,:)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:) = f(:)*j
END DO
END SELECT
END SUBROUTINE ppval1
!===========================================================================
SUBROUTINE ppval1z(sp, x, ppform, left, jder, f)
!
! Compute function and derivatives from the PP representation
! f(x) = \sum_{k=0}^p ppform(k)*(x-t_{left})^k
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE COMPLEX, INTENT(in) :: ppform(:,:)
INTEGER, INTENT(in) :: left, jder
DOUBLE COMPLEX, INTENT(out) :: f(:)
DOUBLE PRECISION :: h, fact
INTEGER :: j, order
!
order = sp%order ! Polynomial degree p + 1
!
h = x-sp%knots(left)
f = (0.0d0,0.0d0)
IF( jder .LT. 0 .OR. jder .GE. order ) RETURN
!
SELECT CASE (jder)
CASE(0) ! function value
DO j=order,1,-1
f(:) = f(:)*h + ppform(j,:)
END DO
CASE(1) ! 1st derivative
DO j=order-1,1,-1
f(:) = f(:)*h + j*ppform(j+1,:)
END DO
CASE default ! 2nd and higher derivatives
fact = order-jder
DO j=order,jder+1,-1
f(:) = f(:)/fact*j*h + ppform(j,:)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:) = f(:)*j
END DO
END SELECT
END SUBROUTINE ppval1z
!===========================================================================
SUBROUTINE ppval2(sp, x, ppform, left, jder, f)
!
! Compute function and derivatives from the PP representation
! f(x) = \sum_{k=0}^p ppform(k)*(x-t_{left})^k
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x, ppform(:,:,:)
INTEGER, INTENT(in) :: left, jder
DOUBLE PRECISION, INTENT(out) :: f(:,:)
DOUBLE PRECISION :: h, fact
INTEGER :: j, order
!
order = sp%order ! Polynomial degree p + 1
!
h = x-sp%knots(left)
f = 0.0d0
IF( jder .LT. 0 .OR. jder .GE. order ) RETURN
!
SELECT CASE (jder)
CASE(0) ! function value
DO j=order,1,-1
f(:,:) = f(:,:)*h + ppform(j,:,:)
END DO
CASE(1) ! 1st derivative
DO j=order-1,1,-1
f(:,:) = f(:,:)*h + j*ppform(j+1,:,:)
END DO
CASE default ! 2nd and higher derivatives
fact = order-jder
DO j=order,jder+1,-1
f(:,:) = f(:,:)/fact*j*h + ppform(j,:,:)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:,:) = f(:,:)*j
END DO
END SELECT
END SUBROUTINE ppval2
!===========================================================================
SUBROUTINE ppval2z(sp, x, ppform, left, jder, f)
!
! Compute function and derivatives from the PP representation
! f(x) = \sum_{k=0}^p ppform(k)*(x-t_{left})^k
!
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE COMPLEX, INTENT(in) :: ppform(:,:,:)
INTEGER, INTENT(in) :: left, jder
DOUBLE COMPLEX, INTENT(out) :: f(:,:)
DOUBLE PRECISION :: h, fact
INTEGER :: j, order
!
order = sp%order ! Polynomial degree p + 1
!
h = x-sp%knots(left)
f = (0.0d0,0.0d0)
IF( jder .LT. 0 .OR. jder .GE. order ) RETURN
!
SELECT CASE (jder)
CASE(0) ! function value
DO j=order,1,-1
f(:,:) = f(:,:)*h + ppform(j,:,:)
END DO
CASE(1) ! 1st derivative
DO j=order-1,1,-1
f(:,:) = f(:,:)*h + j*ppform(j+1,:,:)
END DO
CASE default ! 2nd and higher derivatives
fact = order-jder
DO j=order,jder+1,-1
f(:,:) = f(:,:)/fact*j*h + ppform(j,:,:)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:,:) = f(:,:)*j
END DO
END SELECT
END SUBROUTINE ppval2z
!===========================================================================
SUBROUTINE locintv0_old(sp, x, left)
!
! Locate the interval containing x
! Should be in [0, nints-1]
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: left
DOUBLE PRECISION :: hinv
INTEGER :: nints
!
nints = sp%nints
!
! Case of equidistant mesh
IF( sp%nlequid) THEN
hinv = sp%hinv
left = MAX(0,MIN(FLOOR((x-sp%knots(0))*hinv),nints-1))
RETURN
END IF
!
! Non-equistant mesh
left = sp%left
DO
IF( left .EQ. nints ) THEN
left = nints-1
EXIT
END IF
IF( left .LT. 0 ) THEN
left = 0
EXIT
END IF
IF( x .LT. sp%knots(left+1) ) THEN
IF( x .GE. sp%knots(left) ) THEN
EXIT
ELSE
left = left-1
END IF
ELSE
left = left+1
END IF
END DO
IF(left .GT. 0 .AND. left .LT. nints) THEN
sp%left = left
END IF
END SUBROUTINE locintv0_old
!
!===========================================================================
SUBROUTINE locintv1_old(sp, x, left)
!
! Locate the intervals left(:) containing x(:)
! Should be in [0, nints-1]
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: x(:)
INTEGER, INTENT(out) :: left(:)
DOUBLE PRECISION :: hinv
INTEGER :: nints, i
!
! Case of equidistant mesh
nints = sp%nints
IF( sp%nlequid) THEN
hinv = sp%hinv
left(:) = MAX(0,MIN(FLOOR((x(:)-sp%knots(0))*hinv),nints-1))
RETURN
END IF
!
! Non-equistant mesh
DO i=1,SIZE(x)
CALL locintv0_old(sp, x(i), left(i))
END DO
END SUBROUTINE locintv1_old
!
!===========================================================================
SUBROUTINE locintv0(sp, x, left)
!
! Locate the interval containing x
! Should be in [0, nints-1]
!
- TYPE(spline1d) :: sp
+ TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: left
DOUBLE PRECISION :: hinv
INTEGER :: l, nf, nints
!
nints = sp%nints
!
! Case of equidistant mesh
IF( sp%nlequid) THEN
hinv = sp%hinv
left = MAX(0,MIN(FLOOR((x-sp%knots(0))*hinv),nints-1))
RETURN
END IF
!
! Non-equistant mesh
hinv = sp%hinv
nf = SIZE(sp%fmap) - 1
l = MAX(0,MIN(FLOOR((x-sp%knots(0))*hinv),nf-1)) ! left on fine mesh
left = sp%fmap(l)
IF( x.GE.sp%knots(left+1) ) left = MIN(left+1,nints-1)
END SUBROUTINE locintv0
!
!===========================================================================
SUBROUTINE locintv1(sp, x, left)
!
! Locate the intervals left(:) containing x(:)
! Should be in [0, nints-1]
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: x(:)
INTEGER, INTENT(out) :: left(:)
INTEGER :: l(SIZE(x))
DOUBLE PRECISION :: hinv
INTEGER :: npt, nf, nints, i
!
npt = SIZE(x)
!
! Case of equidistant mesh
nints = sp%nints
IF( sp%nlequid) THEN
hinv = sp%hinv
left(1:npt) = MAX(0,MIN(FLOOR((x(1:npt)-sp%knots(0))*hinv),nints-1))
RETURN
END IF
!
! Non-equistant mesh
hinv = sp%hinv
nf = SIZE(sp%fmap) - 1
l(:) = MAX(0,MIN(FLOOR((x(:)-sp%knots(0))*hinv),nf-1)) ! left on fine mesh
left(1:npt) = sp%fmap(l(1:npt))
WHERE( x.GE.sp%knots(left+1) ) left = MIN(left+1,nints-1)
END SUBROUTINE locintv1
!
!===========================================================================
SUBROUTINE destroy_sp1d(sp)
!
! Clean up 1d spline object
!
TYPE(spline1d) :: sp
!
IF( ASSOCIATED(sp%knots) ) DEALLOCATE (sp%knots)
IF( ASSOCIATED(sp%val0) ) DEALLOCATE (sp%val0)
IF( ASSOCIATED(sp%valc) ) DEALLOCATE (sp%valc)
IF( ASSOCIATED(sp%gausx) ) DEALLOCATE (sp%gausx)
IF( ASSOCIATED(sp%gausw) ) DEALLOCATE (sp%gausw)
IF( ASSOCIATED(sp%intspl) ) DEALLOCATE (sp%intspl)
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE (sp%ppform)
IF( ASSOCIATED(sp%ppformz) ) DEALLOCATE (sp%ppformz)
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
IF( ASSOCIATED(sp%bcoefsc) ) DEALLOCATE(sp%bcoefsc)
IF( ASSOCIATED(sp%fmap) ) DEALLOCATE(sp%fmap)
!
CALL destroy(sp%mat)
CALL destroy(sp%matp)
CALL destroy_dftmap(sp%dft)
END SUBROUTINE destroy_sp1d
!
!===========================================================================
SUBROUTINE destroy_dftmap(m)
!
! Clean up DFTMAP
!
TYPE(dftmap) :: m
!
IF(ASSOCIATED(m%coefs)) DEALLOCATE(m%coefs)
END SUBROUTINE destroy_dftmap
!===========================================================================
SUBROUTINE destroy_sp2d(sp)
!
! Clean up 2d spline object
!
TYPE(spline2d) :: sp
!
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
IF( ASSOCIATED(sp%ppformz) ) DEALLOCATE(sp%ppformz)
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
IF( ASSOCIATED(sp%bcoefsc) ) DEALLOCATE(sp%bcoefsc)
CALL destroy_sp1d(sp%sp1)
CALL destroy_sp1d(sp%sp2)
END SUBROUTINE destroy_sp2d
!===========================================================================
SUBROUTINE destroy_sp2d1d(sp)
!
! Clean up 2d1d spline object
!
TYPE(spline2d1d) :: sp
!
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
IF( ASSOCIATED(sp%ppformz) ) DEALLOCATE (sp%ppformz)
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
IF( ASSOCIATED(sp%bcoefsc) ) DEALLOCATE(sp%bcoefsc)
CALL destroy_sp2d(sp%sp12)
CALL destroy_sp1d(sp%sp3)
END SUBROUTINE destroy_sp2d1d
!
!===========================================================================
SUBROUTINE calc_integ0(sp, finteg)
!
! Compute integral of splines
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(out) :: finteg(0:)
DOUBLE PRECISION, ALLOCATABLE :: xg(:), wg(:), fun(:,:)
DOUBLE PRECISION :: x1, x2
INTEGER :: dim, nx, nidbas, ng, i, ig, j, jj, left
!
CALL get_dim(sp, dim, nx, nidbas)
ng = MAX(2, (nidbas+2)/2)
ALLOCATE(xg(ng), wg(ng), fun(0:nidbas,1))
fun = 0.0d0
finteg = 0.0d0
DO i=1,nx ! Loop thru the intervals
left = i
x1 = sp%knots(i-1)
x2 = sp%knots(i)
CALL gauleg(x1, x2, xg, wg, ng)
DO ig=1,ng ! Loop thru Gauss points
CALL basfun(xg(ig), sp, fun, i)
left = i-1
DO j=0,nidbas ! Loop thru the splines [left:left+nidbas]
jj = left+j ! in this interval
IF( sp%period ) jj = MODULO(left+j, nx)
finteg(jj) = finteg(jj) + wg(ig)*fun(j,1)
END DO
END DO
END DO
DEALLOCATE(xg, wg, fun)
END SUBROUTINE calc_integ0
!===========================================================================
SUBROUTINE calc_integn(sp, finteg)
!
! Compute integrals = Int( x^a \Lambda_j(x) ), a=0,1,... n
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(out) :: finteg(0:,0:)
DOUBLE PRECISION, ALLOCATABLE :: xg(:), wg(:), fun(:,:)
DOUBLE PRECISION :: x1, x2, xpow
INTEGER :: dim, nx, nidbas, ng, i, ig, j, k, jj, left
INTEGER :: nord
!
nord = SIZE(finteg,2)-1
CALL get_dim(sp, dim, nx, nidbas)
ng = MAX(2, (nidbas+nord+2)/2)
ALLOCATE(xg(ng), wg(ng), fun(0:nidbas,1))
fun = 0.0d0
finteg = 0.0d0
DO i=1,nx ! Loop thru the intervals
left = i
x1 = sp%knots(i-1)
x2 = sp%knots(i)
CALL gauleg(x1, x2, xg, wg, ng)
DO ig=1,ng ! Loop thru Gauss points
CALL basfun(xg(ig), sp, fun, i)
left=i-1
DO j=0,nidbas ! Loop thru the splines [left:left+nidbas]
jj = left+j ! in this interval
IF( sp%period ) jj = MODULO(left+j, nx)
xpow = wg(ig)*fun(j,1)
DO k=0,nord
finteg(jj,k) = finteg(jj,k) + xpow
xpow = xpow*xg(ig)
END DO
END DO
END DO
END DO
DEALLOCATE(xg, wg, fun)
END SUBROUTINE calc_integn
!===========================================================================
!
DOUBLE PRECISION FUNCTION fintg1(sp, c)
!
! Integral of 1d function from its spline coefs c.
!
TYPE(spline1d) :: sp
DOUBLE PRECISION, INTENT(in) :: c(0:)
INTEGER :: dim
dim = sp%dim
fintg1 = DOT_PRODUCT(sp%intspl(0:dim-1), c(0:dim-1))
END FUNCTION fintg1
!===========================================================================
DOUBLE PRECISION FUNCTION fintg2(sp, c)
!
! Integral of 2d function from its spline coefs c.
!
TYPE(spline2d) :: sp
DOUBLE PRECISION, INTENT(in) :: c(0:,0:)
INTEGER :: dim1, dim2, i, j
dim1 = sp%sp1%dim
dim2 = sp%sp2%dim
fintg2 = 0.0d0
DO j=0,dim2-1
DO i=0,dim1-1
fintg2 = fintg2 + c(i,j)*sp%sp1%intspl(i)*sp%sp2%intspl(j)
END DO
END DO
END FUNCTION fintg2
!===========================================================================
SUBROUTINE gridval2d_2d(sp, xp, yp, fp, jder, c, ppform)
!
! Compute values or jder-th dervivative of f(x,y) from ppform
! of spline sp. Recompute the ppform if the optional spline
! coefficients c are given.
!
! F(I,J) = F(X(I), Y(J))
!
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp
DOUBLE PRECISION, DIMENSION(:,:), INTENT(out) :: fp
INTEGER, INTENT(in) :: jder(2)
DOUBLE PRECISION, DIMENSION(:,:), OPTIONAL, INTENT(in) :: c
DOUBLE PRECISION, DIMENSION(:,:,:,:), OPTIONAL :: ppform
!
INTEGER :: d1, d2, k1, k2, n1, n2
DOUBLE PRECISION, ALLOCATABLE :: work(:,:,:), temp(:)
DOUBLE PRECISION, ALLOCATABLE :: funx(:,:), funy(:,:)
DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp))
INTEGER :: leftx(SIZE(xp)), lefty(SIZE(yp))
INTEGER :: i, j, k, ii, jj
LOGICAL :: nlppform
!
d1 = sp%sp1%dim
d2 = sp%sp2%dim
k1 = sp%sp1%order
k2 = sp%sp2%order
n1 = sp%sp1%nints
n2 = sp%sp2%nints
nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform
!
! Compute PPFORMM/BCOEFS if spline coefs are passed
!
IF( PRESENT(c) ) THEN
IF( nlppform ) THEN
ALLOCATE(work(d2,k1,n1))
CALL topp1(sp%sp1, c , work)
IF(PRESENT(ppform)) THEN
CALL topp2(sp%sp2, work, ppform)
ELSE
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
ALLOCATE(sp%ppform(k1,n1,k2,n2))
CALL topp2(sp%sp2, work, sp%ppform)
END IF
DEALLOCATE(work)
ELSE
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
ALLOCATE(sp%bcoefs(SIZE(c,1),SIZE(c,2)))
sp%bcoefs = c
END IF
END IF
!
! Applly periodicity if required
!
IF( sp%sp1%period ) THEN ! ** Applly periodicity **
x = sp%sp1%knots(0) + MODULO(xp-sp%sp1%knots(0), sp%sp1%lperiod)
ELSE
x = xp
END IF
IF( sp%sp2%period ) THEN ! ** Applly periodicity **
y = sp%sp2%knots(0) + MODULO(yp-sp%sp2%knots(0), sp%sp2%lperiod)
ELSE
y = yp
END IF
!
! Locate interval containing (x,y)
!
CALL locintv(sp%sp1, x, leftx)
CALL locintv(sp%sp2, y, lefty)
!
! Compute function/derivatives
!
IF( nlppform ) THEN ! using PP form
ALLOCATE(temp(k2))
DO j=1,SIZE(y)
DO i=1,SIZE(x)
IF(PRESENT(ppform)) THEN
CALL ppval(sp%sp1, x(i), ppform(:,leftx(i)+1,:,lefty(j)+1),&
& leftx(i), jder(1), temp)
ELSE
CALL ppval(sp%sp1, x(i), sp%ppform(:,leftx(i)+1,:,lefty(j)+1),&
& leftx(i), jder(1), temp)
END IF
CALL ppval(sp%sp2, y(j), temp, lefty(j), jder(2), fp(i,j))
END DO
END DO
DEALLOCATE(temp)
ELSE ! using spline expansion
ALLOCATE(funx(1:k1,0:jder(1)), funy(1:k2,0:jder(2)))
fp = 0.0d0
DO j=1,SIZE(y)
CALL basfun(y(j), sp%sp2, funy, lefty(j)+1)
DO i=1,SIZE(x)
CALL basfun(x(i), sp%sp1, funx, leftx(i)+1)
DO jj=1,k2
DO ii=1,k1
fp(i,j) = fp(i,j) + sp%bcoefs(leftx(i)+ii,lefty(j)+jj) * &
& funx(ii,jder(1))*funy(jj,jder(2))
END DO
END DO
END DO
END DO
DEALLOCATE(funx, funy)
END IF
END SUBROUTINE gridval2d_2d
!===========================================================================
SUBROUTINE gridval2d1d_3d(sp, xp, yp, zp, fp, jder, c, ppform)
!
! Compute values or jder-th dervivative of f(x,y,z) from spline
! coefficients (nlppform=.false.) or ppform (nlppform=.true.)
!
! F(I,J,K) = F(X(I), Y(J), Z(K))
!
TYPE(spline2d1d), TARGET :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp, zp
DOUBLE PRECISION, DIMENSION(:,:,:), INTENT(out) :: fp
INTEGER, INTENT(in) :: jder(3)
DOUBLE PRECISION, DIMENSION(:,:,:), OPTIONAL, INTENT(in) :: c
DOUBLE PRECISION, DIMENSION(:,:,:,:,:,:), OPTIONAL :: ppform
!
TYPE(spline2d), POINTER :: sp2
DOUBLE PRECISION, ALLOCATABLE :: work1(:,:,:,:), work2(:,:,:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: temp1(:,:), temp2(:)
DOUBLE PRECISION, ALLOCATABLE :: funx(:,:), funy(:,:), funz(:,:)
DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp)), z(SIZE(zp))
INTEGER :: leftx(SIZE(xp)), lefty(SIZE(yp)), leftz(SIZE(zp))
INTEGER :: d1, d2, d3, k1, k2, k3, n1, n2, n3
INTEGER :: ipx, ipy, ipz, k, ii, jj, kk
LOGICAL :: nlppform
!--------------------------------------------------------------------------------
! 1. Prologue
!
sp2 => sp%sp12
d1 = sp2%sp1%dim
d2 = sp2%sp2%dim
d3 = sp%sp3%dim
k1 = sp2%sp1%order
k2 = sp2%sp2%order
k3 = sp%sp3%order
n1 = sp2%sp1%nints
n2 = sp2%sp2%nints
n3 = sp%sp3%nints
nlppform = sp2%sp1%nlppform .OR. sp2%sp2%nlppform .OR. sp%sp3%nlppform
!
! Applly periodicity if required
IF( sp2%sp1%period ) THEN
x = sp2%sp1%knots(0) + MODULO(xp-sp2%sp1%knots(0), sp2%sp1%lperiod)
ELSE
x = xp
END IF
IF( sp2%sp2%period ) THEN
y = sp2%sp2%knots(0) + MODULO(yp-sp2%sp2%knots(0), sp2%sp2%lperiod)
ELSE
y = yp
END IF
IF( sp%sp3%period ) THEN
z = sp%sp3%knots(0) + MODULO(zp-sp%sp3%knots(0), sp%sp3%lperiod)
ELSE
z = zp
END IF
!
! Locate interval containing (x,y,z)
CALL locintv(sp2%sp1, x, leftx)
CALL locintv(sp2%sp2, y, lefty)
CALL locintv(sp%sp3, z, leftz)
!--------------------------------------------------------------------------------
! 2. Using PPFORM
!
IF( nlppform ) THEN
!
! Compute PPFORM from BCOEF
IF( PRESENT(c) ) THEN
ALLOCATE(work2(d3,k1,n1,k2,n2))
ALLOCATE(work1(d2,d3,k1,n1))
CALL topp2(sp2%sp1, c, work1)
CALL topp3(sp2%sp2, work1, work2)
DEALLOCATE(work1)
IF( PRESENT(ppform) )THEN
CALL topp4(sp%sp3, work2, ppform)
ELSE
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
ALLOCATE(sp%ppform(k1,n1,k2,n2,k3,n3))
CALL topp4(sp%sp3, work2, sp%ppform)
END IF
DEALLOCATE(work2)
END IF
!
! Compute function/derivatives
ALLOCATE(temp1(k2,k3))
ALLOCATE(temp2(k3))
DO ipz=1,SIZE(z)
DO ipy=1,SIZE(y)
DO ipx=1,SIZE(x)
IF(PRESENT(ppform)) THEN
CALL ppval(sp2%sp1, x(ipx), &
& ppform(:,leftx(ipx)+1,:,lefty(ipy)+1,:,leftz(ipz)+1),&
& leftx(ipx), jder(1), temp1)
ELSE
CALL ppval(sp2%sp1, x(ipx), &
& sp%ppform(:,leftx(ipx)+1,:,lefty(ipy)+1,:,leftz(ipz)+1),&
& leftx(ipx), jder(1), temp1)
END IF
CALL ppval(sp2%sp2, y(ipy), temp1, lefty(ipy), jder(2), &
& temp2)
CALL ppval(sp%sp3, z(ipz), temp2, leftz(ipz), jder(3), &
& fp(ipx,ipy,ipz))
END DO
END DO
END DO
DEALLOCATE(temp1)
DEALLOCATE(temp2)
!--------------------------------------------------------------------------------
! 3. Using spline expansion
!
ELSE
IF( PRESENT(c) ) THEN
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
ALLOCATE(sp%bcoefs(SIZE(c,1),SIZE(c,2),SIZE(c,3)))
sp%bcoefs = c
END IF
!
! Compute function/derivatives
ALLOCATE(funx(1:k1,0:jder(1)))
ALLOCATE(funy(1:k2,0:jder(2)))
ALLOCATE(funz(1:k3,0:jder(3)))
fp = 0.0d0
DO ipz=1,SIZE(z)
CALL basfun(z(ipz), sp%sp3, funz, leftz(ipz)+1)
DO ipy=1,SIZE(y)
CALL basfun(y(ipy), sp2%sp2, funy, lefty(ipy)+1)
DO ipx=1,SIZE(x)
CALL basfun(x(ipx), sp2%sp1, funx, leftx(ipx)+1)
DO kk=1,k3
DO jj=1,k2
DO ii=1,k1
fp(ipx,ipy,ipz) = fp(ipx,ipy,ipz) + &
& sp%bcoefs(leftx(ipx)+ii,lefty(ipy)+jj,leftz(ipz)+kk) * &
& funx(ii,jder(1)) * funy(jj,jder(2)) * funz(kk,jder(3))
END DO
END DO
END DO
END DO
END DO
END DO
DEALLOCATE(funx, funy, funz)
END IF
END SUBROUTINE gridval2d1d_3d
!===========================================================================
SUBROUTINE gridval2d1d_1d(sp, xp, yp, zp, fp, jder, c, ppform)
!
! Compute values or jder-th dervivative of f(x,y,z) from spline
! coefficients (nlppform=.false.) or ppform (nlppform=.true.)
!
! F(I) = F(X(I),Y(I),Z(I))
!
TYPE(spline2d1d), TARGET :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp, zp
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: fp
INTEGER, INTENT(in) :: jder(3)
DOUBLE PRECISION, DIMENSION(:,:,:), OPTIONAL, INTENT(in) :: c
DOUBLE PRECISION, DIMENSION(:,:,:,:,:,:), OPTIONAL :: ppform
!
TYPE(spline2d), POINTER :: sp2
DOUBLE PRECISION, ALLOCATABLE :: work1(:,:,:,:), work2(:,:,:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: temp1(:,:), temp2(:)
DOUBLE PRECISION, ALLOCATABLE :: funx(:,:), funy(:,:), funz(:,:)
DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp)), z(SIZE(zp))
INTEGER :: leftx(SIZE(xp)), lefty(SIZE(yp)), leftz(SIZE(zp))
INTEGER :: d1, d2, d3, k1, k2, k3, n1, n2, n3
INTEGER :: np, ip, ii, jj, kk
LOGICAL :: nlppform
!--------------------------------------------------------------------------------
! 1. Prologue
!
sp2 => sp%sp12
d1 = sp2%sp1%dim
d2 = sp2%sp2%dim
d3 = sp%sp3%dim
k1 = sp2%sp1%order
k2 = sp2%sp2%order
k3 = sp%sp3%order
n1 = sp2%sp1%nints
n2 = sp2%sp2%nints
n3 = sp%sp3%nints
np = SIZE(xp)
nlppform = sp2%sp1%nlppform .OR. sp2%sp2%nlppform .OR. sp%sp3%nlppform
!
! Applly periodicity if required
IF( sp2%sp1%period ) THEN
x = sp2%sp1%knots(0) + MODULO(xp-sp2%sp1%knots(0), sp2%sp1%lperiod)
ELSE
x = xp
END IF
IF( sp2%sp2%period ) THEN
y = sp2%sp2%knots(0) + MODULO(yp-sp2%sp2%knots(0), sp2%sp2%lperiod)
ELSE
y = yp
END IF
IF( sp%sp3%period ) THEN
z = sp%sp3%knots(0) + MODULO(zp-sp%sp3%knots(0), sp%sp3%lperiod)
ELSE
z = zp
END IF
!
! Locate interval containing (x,y,z)
CALL locintv(sp2%sp1, x, leftx)
CALL locintv(sp2%sp2, y, lefty)
CALL locintv(sp%sp3, z, leftz)
!--------------------------------------------------------------------------------
! 2. Using PPFORM
!
IF( nlppform ) THEN
!
! Compute PPFORM from BCOEF
IF( PRESENT(c) ) THEN
ALLOCATE(work2(d3,k1,n1,k2,n2))
ALLOCATE(work1(d2,d3,k1,n1))
CALL topp2(sp2%sp1, c, work1)
CALL topp3(sp2%sp2, work1, work2)
DEALLOCATE(work1)
IF( PRESENT(ppform) )THEN
CALL topp4(sp%sp3, work2, ppform)
ELSE
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
ALLOCATE(sp%ppform(k1,n1,k2,n2,k3,n3))
CALL topp4(sp%sp3, work2, sp%ppform)
END IF
DEALLOCATE(work2)
END IF
!
! Compute function/derivatives
ALLOCATE(temp1(k2,k3))
ALLOCATE(temp2(k3))
DO ip=1,np
IF(PRESENT(ppform)) THEN
CALL ppval(sp2%sp1, x(ip), &
& ppform(:,leftx(ip)+1,:,lefty(ip)+1,:,leftz(ip)+1),&
& leftx(ip), jder(1), temp1)
ELSE
CALL ppval(sp2%sp1, x(ip), &
& sp%ppform(:,leftx(ip)+1,:,lefty(ip)+1,:,leftz(ip)+1),&
& leftx(ip), jder(1), temp1)
END IF
CALL ppval(sp2%sp2, y(ip), temp1, lefty(ip), jder(2), temp2)
CALL ppval(sp%sp3, z(ip), temp2, leftz(ip), jder(3), fp(ip))
END DO
DEALLOCATE(temp1)
DEALLOCATE(temp2)
!--------------------------------------------------------------------------------
! 3. Using spline expansion
!
ELSE
IF( PRESENT(c) ) THEN
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
ALLOCATE(sp%bcoefs(SIZE(c,1),SIZE(c,2),SIZE(c,3)))
sp%bcoefs = c
END IF
!
! Compute function/derivatives
ALLOCATE(funx(1:k1,0:jder(1)))
ALLOCATE(funy(1:k2,0:jder(2)))
ALLOCATE(funz(1:k3,0:jder(3)))
fp = 0.0d0
DO ip=1,np
CALL basfun(x(ip), sp2%sp1, funx, leftx(ip)+1)
CALL basfun(y(ip), sp2%sp2, funy, lefty(ip)+1)
CALL basfun(z(ip), sp%sp3, funz, leftz(ip)+1)
DO kk=1,k3
DO jj=1,k2
DO ii=1,k1
fp(ip) = fp(ip) + &
& sp%bcoefs(leftx(ip)+ii,lefty(ip)+jj,leftz(ip)+kk) * &
& funx(ii,jder(1))*funy(jj,jder(2))*funz(kk,jder(3))
END DO
END DO
END DO
END DO
DEALLOCATE(funx, funy, funz)
END IF
END SUBROUTINE gridval2d1d_1d
!===========================================================================
SUBROUTINE gridval2d(sp, xp, yp, fp, jder, c, ppform)
!
! Compute values or jder-th dervivative of f(x,y) from ppform
! of spline sp. Recompute the ppform if the optional spline
! coefficients c are given.
!
! F = F(X, Y)
!
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, INTENT(in) :: xp, yp
DOUBLE PRECISION, INTENT(out) :: fp
INTEGER, INTENT(in) :: jder(2)
DOUBLE PRECISION, DIMENSION(:,:), OPTIONAL, INTENT(in) :: c
DOUBLE PRECISION, DIMENSION(:,:,:,:), OPTIONAL :: ppform
!
INTEGER :: d1, d2, k1, k2, n1, n2
DOUBLE PRECISION, ALLOCATABLE :: work(:,:,:), temp(:)
DOUBLE PRECISION, ALLOCATABLE :: funx(:,:), funy(:,:)
DOUBLE PRECISION :: x, y
INTEGER :: leftx, lefty
INTEGER :: i, j, k, ii, jj
LOGICAL :: nlppform
!
d1 = sp%sp1%dim
d2 = sp%sp2%dim
k1 = sp%sp1%order
k2 = sp%sp2%order
n1 = sp%sp1%nints
n2 = sp%sp2%nints
nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform
!
! Compute PPFORM/BCOEFS if spline coefs are passed
!
IF( PRESENT(c)) THEN
IF( nlppform ) THEN
ALLOCATE(work(d2,k1,n1))
CALL topp1(sp%sp1, c , work)
IF(PRESENT(ppform)) THEN
CALL topp2(sp%sp2, work, ppform)
ELSE
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
ALLOCATE(sp%ppform(k1,n1,k2,n2))
CALL topp2(sp%sp2, work, sp%ppform)
END IF
DEALLOCATE(work)
ELSE
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
ALLOCATE(sp%bcoefs(SIZE(c,1),SIZE(c,2)))
sp%bcoefs = c
END IF
END IF
!
! Applly periodicity if required
!
IF( sp%sp1%period ) THEN ! ** Applly periodicity **
x = sp%sp1%knots(0) + MODULO(xp-sp%sp1%knots(0), sp%sp1%lperiod)
ELSE
x = xp
END IF
IF( sp%sp2%period ) THEN ! ** Applly periodicity **
y = sp%sp2%knots(0) + MODULO(yp-sp%sp2%knots(0), sp%sp2%lperiod)
ELSE
y = yp
END IF
!
! Locate the interval containing x, y
!
CALL locintv(sp%sp1, x, leftx)
CALL locintv(sp%sp2, y, lefty)
!
! Compute function/derivatives
!
IF( nlppform ) THEN ! using PP form
ALLOCATE(temp(k2))
IF(PRESENT(ppform)) THEN
CALL ppval(sp%sp1, x, ppform(:,leftx+1,:,lefty+1),&
& leftx, jder(1), temp)
ELSE
CALL ppval(sp%sp1, x, sp%ppform(:,leftx+1,:,lefty+1),&
& leftx, jder(1), temp)
END IF
CALL ppval(sp%sp2, y, temp, lefty, jder(2), fp)
DEALLOCATE(temp)
ELSE ! using spline expansion
ALLOCATE(funx(1:k1,0:jder(1)), funy(1:k2,0:jder(2)))
fp = 0.0d0
CALL basfun(x, sp%sp1, funx, leftx+1)
CALL basfun(y, sp%sp2, funy, lefty+1)
DO jj=1,k2
DO ii=1,k1
fp = fp + &
& funy(jj,jder(2))*sp%bcoefs(leftx+ii,lefty+jj)* &
& funx(ii,jder(1))
END DO
END DO
DEALLOCATE(funx, funy)
END IF
END SUBROUTINE gridval2d
!===========================================================================
SUBROUTINE gridval2d_1d(sp, xp, yp, fp, jder, c, ppform)
!
! Compute values or jder-th dervivative of f(x,y) from ppform
! of spline sp. Recompute the ppform if the optional spline
! coefficients c are given.
!
! F(I) = F(X(I), Y(I))
!
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: fp
INTEGER, INTENT(in) :: jder(2)
DOUBLE PRECISION, DIMENSION(:,:), OPTIONAL, INTENT(in) :: c
DOUBLE PRECISION, DIMENSION(:,:,:,:), OPTIONAL :: ppform
!
INTEGER :: d1, d2, k1, k2, n1, n2, np
DOUBLE PRECISION, ALLOCATABLE :: work(:,:,:), temp(:)
DOUBLE PRECISION, ALLOCATABLE :: funx(:,:), funy(:,:)
DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp))
INTEGER :: leftx(SIZE(xp)), lefty(SIZE(yp))
INTEGER :: i, j, k, ii, jj
LOGICAL :: nlppform
!
d1 = sp%sp1%dim
d2 = sp%sp2%dim
k1 = sp%sp1%order
k2 = sp%sp2%order
n1 = sp%sp1%nints
n2 = sp%sp2%nints
nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform
!
! Compute PPFORM/BCOEFS if spline coefs are passed
!
IF( PRESENT(c)) THEN
IF( nlppform ) THEN
ALLOCATE(work(d2,k1,n1))
CALL topp1(sp%sp1, c , work)
IF(PRESENT(ppform)) THEN
CALL topp2(sp%sp2, work, ppform)
ELSE
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
ALLOCATE(sp%ppform(k1,n1,k2,n2))
CALL topp2(sp%sp2, work, sp%ppform)
END IF
DEALLOCATE(work)
ELSE
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
ALLOCATE(sp%bcoefs(SIZE(c,1),SIZE(c,2)))
sp%bcoefs = c
END IF
END IF
!
! Applly periodicity if required
!
np = SIZE(xp)
IF( sp%sp1%period ) THEN ! ** Applly periodicity **
x = sp%sp1%knots(0) + MODULO(xp-sp%sp1%knots(0), sp%sp1%lperiod)
ELSE
x = xp
END IF
IF( sp%sp2%period ) THEN ! ** Applly periodicity **
y = sp%sp2%knots(0) + MODULO(yp-sp%sp2%knots(0), sp%sp2%lperiod)
ELSE
y = yp
END IF
!
! Locate the interval containing x, y
!
CALL locintv(sp%sp1, x, leftx)
CALL locintv(sp%sp2, y, lefty)
!
! Compute function/derivatives
!
IF( nlppform ) THEN ! using PP form
ALLOCATE(temp(k2))
DO i=1,np
IF(PRESENT(ppform)) THEN
CALL ppval(sp%sp1, x(i), ppform(:,leftx(i)+1,:,lefty(i)+1),&
& leftx(i), jder(1), temp)
ELSE
CALL ppval(sp%sp1, x(i), sp%ppform(:,leftx(i)+1,:,lefty(i)+1),&
& leftx(i), jder(1), temp)
END IF
CALL ppval(sp%sp2, y(i), temp, lefty(i), jder(2), fp(i))
END DO
DEALLOCATE(temp)
ELSE ! using spline expansion
ALLOCATE(funx(1:k1,0:jder(1)), funy(1:k2,0:jder(2)))
fp = 0.0d0
DO i=1,np
CALL basfun(x(i), sp%sp1, funx, leftx(i)+1)
CALL basfun(y(i), sp%sp2, funy, lefty(i)+1)
DO jj=1,k2
DO ii=1,k1
fp(i) = fp(i) + &
& funy(jj,jder(2))*sp%bcoefs(leftx(i)+ii,lefty(i)+jj)* &
& funx(ii,jder(1))
END DO
END DO
END DO
DEALLOCATE(funx, funy)
END IF
END SUBROUTINE gridval2d_1d
!===========================================================================
SUBROUTINE gridval2dz(sp, xp, yp, fp, jder, c, ppformz)
!
! Compute values or jder-th dervivative of f(x,y) from ppform
! of spline sp. Recompute the ppform if the optional spline
! coefficients c are given.
!
! F = F(X, Y)
!
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, INTENT(in) :: xp, yp
DOUBLE COMPLEX, INTENT(out) :: fp
INTEGER, INTENT(in) :: jder(2)
DOUBLE COMPLEX, DIMENSION(:,:), OPTIONAL, INTENT(in) :: c
DOUBLE COMPLEX, DIMENSION(:,:,:,:), OPTIONAL :: ppformz
!
INTEGER :: d1, d2, k1, k2, n1, n2
DOUBLE COMPLEX, ALLOCATABLE :: work(:,:,:), temp(:)
DOUBLE PRECISION, ALLOCATABLE :: funx(:,:), funy(:,:)
DOUBLE PRECISION :: x, y
INTEGER :: leftx, lefty
INTEGER :: i, j, k, ii, jj
LOGICAL :: nlppform
!
d1 = sp%sp1%dim
d2 = sp%sp2%dim
k1 = sp%sp1%order
k2 = sp%sp2%order
n1 = sp%sp1%nints
n2 = sp%sp2%nints
nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform
!
! Compute PPFORM/BCOEFS if spline coefs are passed
!
IF( PRESENT(c)) THEN
IF( nlppform ) THEN
ALLOCATE(work(d2,k1,n1))
CALL topp1z(sp%sp1, c , work)
IF(PRESENT(ppformz)) THEN
CALL topp2z(sp%sp2, work, ppformz)
ELSE
IF( ASSOCIATED(sp%ppformz) ) DEALLOCATE(sp%ppformz)
ALLOCATE(sp%ppformz(k1,n1,k2,n2))
CALL topp2z(sp%sp2, work, sp%ppformz)
END IF
DEALLOCATE(work)
ELSE
IF( ASSOCIATED(sp%bcoefsc) ) DEALLOCATE(sp%bcoefsc)
ALLOCATE(sp%bcoefsc(SIZE(c,1),SIZE(c,2)))
sp%bcoefsc = c
END IF
END IF
!
! Applly periodicity if required
!
IF( sp%sp1%period ) THEN ! ** Applly periodicity **
x = sp%sp1%knots(0) + MODULO(xp-sp%sp1%knots(0), sp%sp1%lperiod)
ELSE
x = xp
END IF
IF( sp%sp2%period ) THEN ! ** Applly periodicity **
y = sp%sp2%knots(0) + MODULO(yp-sp%sp2%knots(0), sp%sp2%lperiod)
ELSE
y = yp
END IF
!
! Locate the interval containing x, y
!
CALL locintv(sp%sp1, x, leftx)
CALL locintv(sp%sp2, y, lefty)
!
! Compute function/derivatives
!
IF( nlppform ) THEN ! using PP form
ALLOCATE(temp(k2))
IF(PRESENT(ppformz)) THEN
CALL ppval(sp%sp1, x, ppformz(:,leftx+1,:,lefty+1),&
& leftx, jder(1), temp)
ELSE
CALL ppval(sp%sp1, x, sp%ppformz(:,leftx+1,:,lefty+1),&
& leftx, jder(1), temp)
END IF
CALL ppval(sp%sp2, y, temp, lefty, jder(2), fp)
DEALLOCATE(temp)
ELSE ! using spline expansion
ALLOCATE(funx(1:k1,0:jder(1)), funy(1:k2,0:jder(2)))
fp = (0.0d0,0.0d0)
CALL basfun(x, sp%sp1, funx, leftx+1)
CALL basfun(y, sp%sp2, funy, lefty+1)
DO jj=1,k2
DO ii=1,k1
fp = fp + &
& funy(jj,jder(2))*sp%bcoefsc(leftx+ii,lefty+jj)* &
& funx(ii,jder(1))
END DO
END DO
DEALLOCATE(funx, funy)
END IF
END SUBROUTINE gridval2dz
!===========================================================================
SUBROUTINE gridval2d_1dz(sp, xp, yp, fp, jder, c, ppformz)
!
! Compute values or jder-th dervivative of f(x,y) from ppform
! of spline sp. Recompute the ppform if the optional spline
! coefficients c are given.
!
! F(I) = F(X(I), Y(I))
!
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: fp
INTEGER, INTENT(in) :: jder(2)
DOUBLE COMPLEX, DIMENSION(:,:), OPTIONAL, INTENT(in) :: c
DOUBLE COMPLEX, DIMENSION(:,:,:,:), OPTIONAL :: ppformz
!
INTEGER :: d1, d2, k1, k2, n1, n2, np
DOUBLE COMPLEX, ALLOCATABLE :: work(:,:,:), temp(:)
DOUBLE PRECISION, ALLOCATABLE :: funx(:,:), funy(:,:)
DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp))
INTEGER :: leftx(SIZE(xp)), lefty(SIZE(yp))
INTEGER :: i, j, k, ii, jj
LOGICAL :: nlppform
!
d1 = sp%sp1%dim
d2 = sp%sp2%dim
k1 = sp%sp1%order
k2 = sp%sp2%order
n1 = sp%sp1%nints
n2 = sp%sp2%nints
nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform
!
! Compute PPFORM/BCOEFS if spline coefs are passed
!
IF( PRESENT(c)) THEN
IF( nlppform ) THEN
ALLOCATE(work(d2,k1,n1))
CALL topp1z(sp%sp1, c , work)
IF(PRESENT(ppformz)) THEN
CALL topp2z(sp%sp2, work, ppformz)
ELSE
IF( ASSOCIATED(sp%ppformz) ) DEALLOCATE(sp%ppformz)
ALLOCATE(sp%ppformz(k1,n1,k2,n2))
CALL topp2z(sp%sp2, work, sp%ppformz)
END IF
DEALLOCATE(work)
ELSE
IF( ASSOCIATED(sp%bcoefsc) ) DEALLOCATE(sp%bcoefsc)
ALLOCATE(sp%bcoefsc(SIZE(c,1),SIZE(c,2)))
sp%bcoefsc = c
END IF
END IF
!
! Applly periodicity if required
!
np = SIZE(xp)
IF( sp%sp1%period ) THEN ! ** Applly periodicity **
x = sp%sp1%knots(0) + MODULO(xp-sp%sp1%knots(0), sp%sp1%lperiod)
ELSE
x = xp
END IF
IF( sp%sp2%period ) THEN ! ** Applly periodicity **
y = sp%sp2%knots(0) + MODULO(yp-sp%sp2%knots(0), sp%sp2%lperiod)
ELSE
y = yp
END IF
!
! Locate the interval containing x, y
!
CALL locintv(sp%sp1, x, leftx)
CALL locintv(sp%sp2, y, lefty)
!
! Compute function/derivatives
!
IF( nlppform ) THEN ! using PP form
ALLOCATE(temp(k2))
DO i=1,np
IF(PRESENT(ppformz)) THEN
CALL ppval(sp%sp1, x(i), ppformz(:,leftx(i)+1,:,lefty(i)+1),&
& leftx(i), jder(1), temp)
ELSE
CALL ppval(sp%sp1, x(i), sp%ppformz(:,leftx(i)+1,:,lefty(i)+1),&
& leftx(i), jder(1), temp)
END IF
CALL ppval(sp%sp2, y(i), temp, lefty(i), jder(2), fp(i))
END DO
DEALLOCATE(temp)
ELSE ! using spline expansion
ALLOCATE(funx(1:k1,0:jder(1)), funy(1:k2,0:jder(2)))
fp = (0.0d0,0.0d0)
DO i=1,np
CALL basfun(x(i), sp%sp1, funx, leftx(i)+1)
CALL basfun(y(i), sp%sp2, funy, lefty(i)+1)
DO jj=1,k2
DO ii=1,k1
fp(i) = fp(i) + &
& funy(jj,jder(2))*sp%bcoefsc(leftx(i)+ii,lefty(i)+jj)* &
& funx(ii,jder(1))
END DO
END DO
END DO
DEALLOCATE(funx, funy)
END IF
END SUBROUTINE gridval2d_1dz
!===========================================================================
SUBROUTINE gridval2d_2dz(sp, xp, yp, fp, jder, c, ppformz)
!
! Compute values or jder-th dervivative of f(x,y) from ppform
! of spline sp. Recompute the ppform if the optional spline
! coefficients c are given.
!
! F(I,J) = F(X(I), Y(J))
!
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp
DOUBLE COMPLEX, DIMENSION(:,:), INTENT(out) :: fp
INTEGER, INTENT(in) :: jder(2)
DOUBLE COMPLEX, DIMENSION(:,:), OPTIONAL, INTENT(in) :: c
DOUBLE COMPLEX, DIMENSION(:,:,:,:), OPTIONAL :: ppformz
!
INTEGER :: d1, d2, k1, k2, n1, n2
DOUBLE COMPLEX, ALLOCATABLE :: work(:,:,:), temp(:)
DOUBLE PRECISION, ALLOCATABLE :: funx(:,:), funy(:,:)
DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp))
INTEGER :: leftx(SIZE(xp)), lefty(SIZE(yp))
INTEGER :: i, j, k, ii, jj
LOGICAL :: nlppform
!
d1 = sp%sp1%dim
d2 = sp%sp2%dim
k1 = sp%sp1%order
k2 = sp%sp2%order
n1 = sp%sp1%nints
n2 = sp%sp2%nints
nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform
!
! Compute PPFORMM/BCOEFS if spline coefs are passed
!
IF( PRESENT(c) ) THEN
IF( nlppform ) THEN
ALLOCATE(work(d2,k1,n1))
CALL topp1z(sp%sp1, c , work)
IF(PRESENT(ppformz)) THEN
CALL topp2z(sp%sp2, work, ppformz)
ELSE
IF( ASSOCIATED(sp%ppformz) ) DEALLOCATE(sp%ppformz)
ALLOCATE(sp%ppformz(k1,n1,k2,n2))
CALL topp2z(sp%sp2, work, sp%ppformz)
END IF
DEALLOCATE(work)
ELSE
IF( ASSOCIATED(sp%bcoefsc) ) DEALLOCATE(sp%bcoefsc)
ALLOCATE(sp%bcoefsc(SIZE(c,1),SIZE(c,2)))
sp%bcoefsc = c
END IF
END IF
!
! Applly periodicity if required
!
IF( sp%sp1%period ) THEN ! ** Applly periodicity **
x = sp%sp1%knots(0) + MODULO(xp-sp%sp1%knots(0), sp%sp1%lperiod)
ELSE
x = xp
END IF
IF( sp%sp2%period ) THEN ! ** Applly periodicity **
y = sp%sp2%knots(0) + MODULO(yp-sp%sp2%knots(0), sp%sp2%lperiod)
ELSE
y = yp
END IF
!
! Locate interval containing (x,y)
!
CALL locintv(sp%sp1, x, leftx)
CALL locintv(sp%sp2, y, lefty)
!
! Compute function/derivatives
!
IF( nlppform ) THEN ! using PP form
ALLOCATE(temp(k2))
DO j=1,SIZE(y)
DO i=1,SIZE(x)
IF(PRESENT(ppformz)) THEN
CALL ppval(sp%sp1, x(i), ppformz(:,leftx(i)+1,:,lefty(j)+1),&
& leftx(i), jder(1), temp)
ELSE
CALL ppval(sp%sp1, x(i), sp%ppformz(:,leftx(i)+1,:,lefty(j)+1),&
& leftx(i), jder(1), temp)
END IF
CALL ppval(sp%sp2, y(j), temp, lefty(j), jder(2), fp(i,j))
END DO
END DO
DEALLOCATE(temp)
ELSE ! using spline expansion
ALLOCATE(funx(1:k1,0:jder(1)), funy(1:k2,0:jder(2)))
fp = 0.0d0
DO j=1,SIZE(y)
CALL basfun(y(j), sp%sp2, funy, lefty(j)+1)
DO i=1,SIZE(x)
CALL basfun(x(i), sp%sp1, funx, leftx(i)+1)
DO jj=1,k2
DO ii=1,k1
fp(i,j) = fp(i,j) + sp%bcoefsc(leftx(i)+ii,lefty(j)+jj) * &
& funx(ii,jder(1))*funy(jj,jder(2))
END DO
END DO
END DO
END DO
DEALLOCATE(funx, funy)
END IF
END SUBROUTINE gridval2d_2dz
!===========================================================================
SUBROUTINE calc_fftmass(spl, fftmat)
!
! Compute FT of mass matrix for periodic spline on equidistant mesh
!
TYPE(spline1d) :: spl
DOUBLE PRECISION, INTENT(out) :: fftmat(0:)
!
INTEGER :: dim, nx, nidbas, ngauss
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE COMPLEX, ALLOCATABLE :: ft_fun(:,:)
INTEGER :: igauss, intv
!
CALL get_dim(spl, dim, nx, nidbas)
CALL get_gauss(spl, ngauss)
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
ALLOCATE(ft_fun(0:nx-1,1))
!
! Integrate on first interval
intv = 1
CALL get_gauss(spl, ngauss, intv, xgauss, wgauss)
fftmat = 0.0d0
DO igauss=1,ngauss
CALL ft_basfun(xgauss(igauss), spl, ft_fun, intv)
fftmat(:) = fftmat(:) + wgauss(igauss)*ft_fun(:,1)*CONJG(ft_fun(:,1))
END DO
!
DEALLOCATE(ft_fun)
DEALLOCATE(xgauss, wgauss)
END SUBROUTINE calc_fftmass
!===========================================================================
SUBROUTINE calc_fftmass_old(spl, fftmat)
!
! Compute FT of mass matrix for periodic spline on equidistant mesh
!
TYPE(spline1d) :: spl
DOUBLE PRECISION, INTENT(out) :: fftmat(0:)
INTEGER :: dim, nx, nidbas, ngauss
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:), xgauss(:), wgauss(:), arow(:)
INTEGER :: i, j, k, igauss, intv
DOUBLE PRECISION :: pi, arg0, arg
!
CALL get_dim(spl, dim, nx, nidbas)
CALL get_gauss(spl, ngauss)
ALLOCATE(fun(0:nidbas,1)) ! Spline
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
ALLOCATE(arow(0:nidbas))
!
! Assemble the first row of the upper mass matrix
arow = 0.0d0
intv = 1 ! Get splines on Gauss points in first interval
CALL get_gauss(spl, ngauss, intv, xgauss, wgauss)
DO igauss=1,ngauss
CALL basfun(xgauss(igauss), spl, fun, intv)
DO i=0,nidbas
DO j=0,nidbas-i
arow(i)=arow(i)+fun(j,1)*fun(i+j,1)*wgauss(igauss)
END DO
END DO
END DO
!
! Fourier transform
pi = 4.0d0*ATAN(1.0d0)
arg0 = 2.0d0*pi/REAL(nx,8)
DO k=0,nx-1
fftmat(k) = arow(0)
arg = k*arg0
DO i=1,nidbas
fftmat(k) = fftmat(k) + 2.0d0*arow(i)*COS(i*arg)
END DO
END DO
!
DEALLOCATE(arow)
DEALLOCATE(fun)
DEALLOCATE(xgauss, wgauss)
END SUBROUTINE calc_fftmass_old
!===========================================================================
SUBROUTINE CompMassMatrix1(sp1, sp2, a, b, MassMatrix)
! Compute cross mass matrix MassMatrix between splines sp1 and sp2 over
! interval [a, b]
IMPLICIT NONE
TYPE(spline1d), INTENT(IN) :: sp1, sp2
DOUBLE PRECISION, INTENT(IN) :: a, b
DOUBLE PRECISION, DIMENSION(:, :), POINTER :: MassMatrix
INTEGER :: ndim1, n1, nidbas1
INTEGER :: ndim2, n2, nidbas2
INTEGER :: nint, int, ngauss, ig, k1, k2, i1, j2, left1, left2
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: allknots, xg, wg
DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: fun1, fun2
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(sp1, ndim1, n1, nidbas1)
CALL get_dim(sp2, ndim2, n2, nidbas2)
! PRINT "('In CompMassMatrix1')"
! PRINT "('sp1: dim, #intervals, degree', I, I, I)", ndim1, n1, nidbas1
! PRINT "('sp2: dim, #intervals, degree', I, I, I)", ndim2, n2, nidbas2
ALLOCATE(fun1(0:nidbas1, 1)) ! needs only basis functions (no derivatives)
ALLOCATE(fun2(0:nidbas2, 1)) ! needs only basis functions (no derivatives)
IF (sp1%period) THEN
IF (sp2%period) THEN
ALLOCATE(MassMatrix(n1, n2))
ELSE
ALLOCATE(MassMatrix(n1, ndim2))
END IF
ELSE
IF (sp2%period) THEN
ALLOCATE(MassMatrix(ndim1, n2))
ELSE
ALLOCATE(MassMatrix(ndim1, ndim2))
END IF
END IF
!
! Gauss quadature
!
ALLOCATE(allknots(0:n1+n2+3))
CALL sorted_merge(sp1%knots(0:n1), n1+1, sp2%knots(0:n2), n2+1, a, b, allknots, nint)
nint = nint-1
ngauss = CEILING(REAL(nidbas1 + nidbas2 + 1, 8)/2.D0)
ALLOCATE(xg(ngauss), wg(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
MassMatrix = 0.d0
DO int = 1, nint
! Get gauss abscissas and weights for current interval
CALL gauleg(allknots(int-1), allknots(int), xg, wg, ngauss)
DO ig = 1, ngauss
CALL locintv(sp1, xg(ig), left1)
CALL locintv(sp2, xg(ig), left2)
CALL basfun(xg(ig), sp1, fun1, left1+1)
CALL basfun(xg(ig), sp2, fun2, left2+1)
DO k1 = 0, nidbas1
IF (sp1%period) THEN
i1 = modulo(left1+1 + k1 -1, n1) +1
ELSE
i1 = left1+1 + k1
END IF
DO k2 = 0, nidbas2
IF (sp2%period) THEN
j2 = modulo(left2+1 + k2 -1, n2) +1
ELSE
j2 = left2+1 + k2
END IF
MassMatrix(i1, j2) = MassMatrix(i1, j2) + wg(ig)*fun1(k1, 1)*fun2(k2, 1)
END DO
END DO
END DO
END DO
!===========================================================================
! 3.0 Epilogue
!
DEALLOCATE(xg, wg)
DEALLOCATE(fun1, fun2)
DEALLOCATE(allknots)
END SUBROUTINE CompMassMatrix1
!===========================================================================
SUBROUTINE CompMassMatrix_gb(sp1, sp2, a, b, MassMatrix)
! Compute cross mass matrix MassMatrix between splines sp1 and sp2 over
! interval [a, b]
IMPLICIT NONE
TYPE(spline1d), INTENT(IN) :: sp1, sp2
DOUBLE PRECISION, INTENT(IN) :: a, b
TYPE(gbmat) :: MassMatrix
INTEGER :: ndim1, n1, nidbas1
INTEGER :: ndim2, n2, nidbas2
INTEGER :: nint, int, ngauss, ig, k1, k2, i1, j2, left1, left2
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: allknots, xg, wg
DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: fun1, fun2
DOUBLE PRECISION :: val
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(sp1, ndim1, n1, nidbas1)
CALL get_dim(sp2, ndim2, n2, nidbas2)
ALLOCATE(fun1(0:nidbas1, 1)) ! needs only basis functions (no derivatives)
ALLOCATE(fun2(0:nidbas2, 1)) ! needs only basis functions (no derivatives)
!
! Gauss quadature
!
ALLOCATE(allknots(0:n1+n2+3))
CALL sorted_merge(sp1%knots(0:n1), n1+1, sp2%knots(0:n2), n2+1, a, b, allknots, nint)
nint = nint-1
ngauss = CEILING(REAL(nidbas1 + nidbas2 + 1, 8)/2.D0)
ALLOCATE(xg(ngauss), wg(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
DO int = 1, nint
! Get gauss abscissas and weights for current interval
CALL gauleg(allknots(int-1), allknots(int), xg, wg, ngauss)
DO ig = 1, ngauss
CALL locintv(sp1, xg(ig), left1)
CALL locintv(sp2, xg(ig), left2)
CALL basfun(xg(ig), sp1, fun1, left1+1)
CALL basfun(xg(ig), sp2, fun2, left2+1)
DO k1 = 0, nidbas1
IF (sp1%period) THEN
i1 = modulo(left1+1 + k1 -1, n1) +1
ELSE
i1 = left1+1 + k1
END IF
DO k2 = 0, nidbas2
IF (sp2%period) THEN
j2 = modulo(left2+1 + k2 -1, n2) +1
ELSE
j2 = left2+1 + k2
END IF
val = wg(ig)*fun1(k1, 1)*fun2(k2, 1)
CALL updtmat(MassMatrix, i1, j2, val)
END DO
END DO
END DO
END DO
!===========================================================================
! 3.0 Epilogue
!
DEALLOCATE(xg, wg)
DEALLOCATE(fun1, fun2)
DEALLOCATE(allknots)
END SUBROUTINE CompMassMatrix_gb
!===========================================================================
SUBROUTINE CompMassMatrix_zgb(sp1, sp2, a, b, MassMatrix)
! Compute cross mass matrix MassMatrix between splines sp1 and sp2 over
! interval [a, b]
IMPLICIT NONE
TYPE(spline1d), INTENT(IN) :: sp1, sp2
DOUBLE PRECISION, INTENT(IN) :: a, b
TYPE(zgbmat) :: MassMatrix
INTEGER :: ndim1, n1, nidbas1
INTEGER :: ndim2, n2, nidbas2
INTEGER :: nint, int, ngauss, ig, k1, k2, i1, j2, left1, left2
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: allknots, xg, wg
DOUBLE PRECISION, DIMENSION(:, :), ALLOCATABLE :: fun1, fun2
DOUBLE COMPLEX :: val
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(sp1, ndim1, n1, nidbas1)
CALL get_dim(sp2, ndim2, n2, nidbas2)
ALLOCATE(fun1(0:nidbas1, 1)) ! needs only basis functions (no derivatives)
ALLOCATE(fun2(0:nidbas2, 1)) ! needs only basis functions (no derivatives)
!
! Gauss quadature
!
ALLOCATE(allknots(0:n1+n2+3))
CALL sorted_merge(sp1%knots(0:n1), n1+1, sp2%knots(0:n2), n2+1, a, b, allknots, nint)
nint = nint-1
ngauss = CEILING(REAL(nidbas1 + nidbas2 + 1, 8)/2.D0)
ALLOCATE(xg(ngauss), wg(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
DO int = 1, nint
! Get gauss abscissas and weights for current interval
CALL gauleg(allknots(int-1), allknots(int), xg, wg, ngauss)
DO ig = 1, ngauss
CALL locintv(sp1, xg(ig), left1)
CALL locintv(sp2, xg(ig), left2)
CALL basfun(xg(ig), sp1, fun1, left1+1)
CALL basfun(xg(ig), sp2, fun2, left2+1)
DO k1 = 0, nidbas1
IF (sp1%period) THEN
i1 = modulo(left1+1 + k1 -1, n1) +1
ELSE
i1 = left1+1 + k1
END IF
DO k2 = 0, nidbas2
IF (sp2%period) THEN
j2 = modulo(left2+1 + k2 -1, n2) +1
ELSE
j2 = left2+1 + k2
END IF
val = wg(ig)*fun1(k1, 1)*fun2(k2, 1)
CALL updtmat(MassMatrix, i1, j2, val)
END DO
END DO
END DO
END DO
!===========================================================================
! 3.0 Epilogue
!
DEALLOCATE(xg, wg)
DEALLOCATE(fun1, fun2)
DEALLOCATE(allknots)
END SUBROUTINE CompMassMatrix_zgb
!===========================================================================
SUBROUTINE sorted_merge(arr1, n1, arr2, n2, a, b, arrm, nm)
IMPLICIT NONE
! Peforms:
! 1) Merge of arrays arr1 & arr2 including boundary values a & b
! 2) Sorts the merged arrays keeping only values in [a, b]
! 3) Removes duplicates
INTEGER, INTENT(IN) :: n1, n2
INTEGER, INTENT(OUT) :: nm
DOUBLE PRECISION, INTENT(IN) :: a, b
DOUBLE PRECISION, INTENT(IN) :: arr1(n1), arr2(n2)
DOUBLE PRECISION, DIMENSION(*), INTENT(OUT) :: arrm
INTEGER :: i, j
! Merge the two arrays including a & b
nm = n1 + n2 + 2
arrm(1:nm) = (/ a, arr1(1:n1), b, arr2(1:n2) /)
! Sort
CALL sort(arrm, nm)
! Remove duplicates
j = 1
DO i = 2, nm
IF(arrm(i) .GT. arrm(j)) THEN
j = j + 1
arrm(j) = arrm(i)
END IF
END DO
nm = j
! Remove values outside [a, b]
j = 0
DO i = 1, nm
IF((arrm(i) .GE. a) .AND. (arrm(i) .LE. b)) THEN
j = j + 1
arrm(j) = arrm(i)
END IF
END DO
nm = j
END SUBROUTINE sorted_merge
!===========================================================================
SUBROUTINE sort(arr, n)
! Sorts array ARR of length N into ascending numerical order by the
! Shell-Mezgar algorithm.
! See Sec. 8.1 of Numerical Recipes
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
DOUBLE PRECISION, DIMENSION(n), INTENT(INOUT) :: arr
INTEGER :: nsort, is, i, j, l, m
DOUBLE PRECISION :: tmp
DOUBLE PRECISION, PARAMETER :: tiny = 1D-5
nsort = INT(LOG(REAL(n, 8))/LOG(2.D0) + tiny)
m = n
DO is = 1, nsort
m = m/2
DO j = 1, n-m
i = j
DO
l = i+m
IF (arr(l) .LT. arr(i)) THEN
tmp = arr(i)
arr(i) = arr(l)
arr(l) = tmp
i = i-m
IF (i .LT. 1) EXIT
ELSE
EXIT
END IF
END DO
END DO
END DO
END SUBROUTINE sort
!===========================================================================
LOGICAL FUNCTION is_equid(x, dev)
!
! Check whether mesh is euidistant or not
!
DOUBLE PRECISION, INTENT(in) :: x(0:)
DOUBLE PRECISION, INTENT(out), OPTIONAL :: dev
!
DOUBLE PRECISION :: dx(SIZE(x)-1), dxmin, dxmax, dxaver, e
DOUBLE PRECISION, PARAMETER :: tol=1.d-6
INTEGER :: n, i
n=SIZE(x)-1
dx = (/ (x(i)-x(i-1),i=1,n) /)
dxmin = MINVAL(dx)
dxmax = MAXVAL(dx)
dxaver = (x(n)-x(0))/REAL(n,8)
e = (dxmax-dxmin)/dxaver
!!$ e = (dxmax-dxmin)/(SUM(x)/REAL(n+1))
is_equid = e.LT.tol
IF(PRESENT(dev)) dev = e
END FUNCTION is_equid
!===========================================================================
SUBROUTINE create_fine(cmesh, h, fmap)
!
! Create a fine mesh from a coarse mesh and returns its mapping
!
DOUBLE PRECISION, INTENT(in) :: cmesh(0:)
DOUBLE PRECISION, INTENT(out) :: h
INTEGER, POINTER, INTENT(out) :: fmap(:)
!
DOUBLE PRECISION, ALLOCATABLE :: fmesh(:)
DOUBLE PRECISION :: xlen, hmin
INTEGER :: n, nfine, i, ic
!
n = SIZE(cmesh)-1
xlen = cmesh(n)-cmesh(0)
!
! Minimum interval size
hmin = xlen
DO i=1,n
hmin = MIN(hmin, cmesh(i)-cmesh(i-1))
END DO
!
! Create the fine mesh
nfine = CEILING(xlen/hmin)
h = xlen / REAL(nfine,8)
ALLOCATE(fmap(0:nfine))
ALLOCATE(fmesh(0:nfine))
fmesh = cmesh(0) + (/ (i*h, i=0,nfine) /)
fmesh(nfine) = cmesh(n)
!
! Map fine to coarse mesh
ic = 0
fmap(0) = ic
DO i=1,nfine-1
DO
IF(fmesh(i).GE.cmesh(ic+1)) THEN
ic = ic+1
ELSE
EXIT
END IF
END DO
fmap(i) = ic
END DO
fmap(nfine) = n-1
!
DEALLOCATE(fmesh)
END SUBROUTINE create_fine
!===========================================================================
SUBROUTINE getgradr(sp, xp, yp, f00, f10, f01)
!
! Compute the function f00 and its derivatives
! f10 = d/dx f
! f01 = d/dy f
! assuming that its PPFORM/BCOEFSC was already computed!
!
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp
DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: f00, f10, f01
!
INTEGER :: np
DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp))
INTEGER :: leftx(SIZE(xp)), lefty(SIZE(yp))
INTEGER :: i, ip, ii, jj, nidbas(2)
DOUBLE PRECISION :: temp0(SIZE(xp),sp%sp2%order), temp1(SIZE(xp),sp%sp2%order)
DOUBLE PRECISION, ALLOCATABLE, SAVE :: funx(:,:), funy(:,:)
DOUBLE PRECISION, ALLOCATABLE, SAVE :: ftemp0(:), ftemp1(:)
LOGICAL :: nlppform
!
! Apply periodicity if required
!
np = SIZE(xp)
nidbas(1) = sp%sp1%order-1
nidbas(2) = sp%sp2%order-1
nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform
!
! Locate the interval containing x, y
!
CALL locintv(sp%sp1, xp, leftx)
CALL locintv(sp%sp2, yp, lefty)
x(:) = xp(:) - sp%sp1%knots(leftx(:))
y(:) = yp(:) - sp%sp2%knots(lefty(:))
!
! Compute function/derivatives
!
IF(nlppform) THEN
!
! Using PPFORM
!----------
DO i=1,np
CALL my_ppval1(nidbas(1), x(i), sp%ppform(:,leftx(i)+1,:,lefty(i)+1), &
& temp0(i,:), temp1(i,:))
END DO
!
CALL my_ppval0(nidbas(2), y, temp0, 0, f00)
CALL my_ppval0(nidbas(2), y, temp0, 1, f01)
CALL my_ppval0(nidbas(2), y, temp1, 0, f10)
ELSE
!
! Using spline expansion with sp%bcoefsc
!----------
IF(.NOT.ALLOCATED(funx)) THEN
ALLOCATE(funx(0:nidbas(1),0:1)) ! Spline and its first derivative
ALLOCATE(funy(0:nidbas(2),0:1))
ALLOCATE(ftemp0(0:nidbas(1)))
ALLOCATE(ftemp1(0:nidbas(1)))
END IF
!
DO ip=1,np
CALL my_splines(nidbas(1), x(ip), sp%sp1%val0(:,:,leftx(ip)+1), funx)
CALL my_splines(nidbas(2), y(ip), sp%sp2%val0(:,:,lefty(ip)+1), funy)
DO ii=0,nidbas(1)
ftemp0(ii) = (0.d0,0.d0)
ftemp1(ii) = (0.d0,0.d0)
DO jj=0,nidbas(2)
ftemp0(ii) = ftemp0(ii) + sp%bcoefs(leftx(ip)+1+ii,lefty(ip)+1+jj)*funy(jj,0)
ftemp1(ii) = ftemp1(ii) + sp%bcoefs(leftx(ip)+1+ii,lefty(ip)+1+jj)*funy(jj,1)
END DO
END DO
f00(ip) = SUM(funx(:,0)*ftemp0(:))
f01(ip) = SUM(funx(:,0)*ftemp1(:))
f10(ip) = SUM(funx(:,1)*ftemp0(:))
END DO
!-----------
END IF
CONTAINS
!+++
SUBROUTINE my_ppval0(p, x, ppform, jder, f)
!
! Compute function and derivatives from the PP representation
! for many points x(:)
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION, INTENT(in) :: ppform(:,:)
INTEGER, INTENT(in) :: jder
DOUBLE PRECISION, INTENT(out) :: f(:)
DOUBLE PRECISION :: fact
INTEGER :: j
SELECT CASE (jder)
CASE(0) ! function value
SELECT CASE(p)
CASE(1)
f(:) = ppform(:,1) + x(:)*ppform(:,2)
CASE(2)
f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*ppform(:,3))
!!$ CASE(3)
!!$ f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*(ppform(:,3)+x(:)*ppform(:,4)))
CASE(3:)
f(:) = ppform(:,p+1)
DO j=p,1,-1
f(:) = f(:)*x(:) + ppform(:,j)
END DO
END SELECT
CASE(1) ! 1st derivative
SELECT CASE(p)
CASE(1)
f(:) = ppform(:,2)
CASE(2)
f(:) = ppform(:,2) + x(:)*2.d0*ppform(:,3)
!!$ CASE(3)
!!$ f(:) = ppform(:,2) + x(:)*(2.d0*ppform(:,3)+x(:)*3.0d0*ppform(:,4))
CASE(3:)
f(:) = p*ppform(:,p+1)
DO j=p-1,1,-1
f(:) = f(:)*x(:) + j*ppform(:,j+1)
END DO
END SELECT
CASE default ! 2nd and higher derivatives
f(:) = ppform(:,p+1)
fact = p-jder
DO j=p,jder+1,-1
f(:) = f(:)/fact*j*x(:) + ppform(:,j)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:) = f(:)*j
END DO
END SELECT
END SUBROUTINE my_ppval0
!+++
SUBROUTINE my_ppval1(p, x, ppform, f0, f1)
!
! Compute function and first derivative from the PP representation
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE PRECISION, INTENT(in) :: ppform(:,:)
DOUBLE PRECISION, INTENT(out) :: f0(:)
DOUBLE PRECISION, INTENT(out) :: f1(:)
DOUBLE PRECISION :: fact
INTEGER :: j
SELECT CASE(p)
CASE(1)
f0(:) = ppform(1,:) + x*ppform(2,:)
f1(:) = ppform(2,:)
CASE(2)
f0(:) = ppform(1,:) + x*(ppform(2,:)+x*ppform(3,:))
f1(:) = ppform(2,:) + x*2.d0*ppform(3,:)
CASE(3)
f0(:) = ppform(1,:) + x*(ppform(2,:)+x*(ppform(3,:)+x*ppform(4,:)))
f1(:) = ppform(2,:) + x*(2.d0*ppform(3,:)+x*3.0d0*ppform(4,:))
CASE(4:)
f0 = ppform(p+1,:)
f1 = f0
DO j=p,2,-1
f0(:) = ppform(j,:) + x*f0(:)
f1(:) = f0(:) + x*f1(:)
END DO
f0(:) = ppform(1,:) + x*f0(:)
END SELECT
END SUBROUTINE my_ppval1
!+++
SUBROUTINE my_splines(p, x, ppform, f)
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE PRECISION, INTENT(in) :: ppform(0:p,0:p)
DOUBLE PRECISION, INTENT(out) :: f(0:p,0:1)
INTEGER :: i
DOUBLE PRECISION :: powerx(0:p)
SELECT CASE(p)
CASE(1)
f(0,0) = ppform(0,0) + x*ppform(1,0)
f(0,1) = ppform(1,0)
f(1,0) = 1.0-f(0,0)
f(1,1) = -f(0,1)
CASE(2)
f(0,0) = ppform(0,0) + x*(ppform(1,0)+x*ppform(2,0))
f(0,1) = ppform(1,0) + 2.d0*x*ppform(2,0)
f(1,0) = ppform(0,1) + x*(ppform(1,1)+x*ppform(2,1))
f(1,1) = ppform(1,1) + 2.d0*x*ppform(2,1)
f(2,0) = 1.0 - f(0,0) - f(1,0)
f(2,1) = - f(0,1) - f(1,1)
CASE(3)
f(0,0) = ppform(0,0) + x*(ppform(1,0)+x*(ppform(2,0)+x*ppform(3,0)))
f(0,1) = ppform(1,0) + x*(2.d0*ppform(2,0)+3.d0*x*ppform(3,0))
f(1,0) = ppform(0,1) + x*(ppform(1,1)+x*(ppform(2,1)+x*ppform(3,1)))
f(1,1) = ppform(1,1) + x*(2.d0*ppform(2,1)+3.d0*x*ppform(3,1))
f(2,0) = ppform(0,2) + x*(ppform(1,2)+x*(ppform(2,2)+x*ppform(3,2)))
f(2,1) = ppform(1,2) + x*(2.d0*ppform(2,2)+3.d0*x*ppform(3,2))
f(3,0) = 1.0 - f(0,0) - f(1,0) - f(2,0)
f(3,1) = - f(0,1) - f(1,1) - f(2,1)
CASE(4:)
powerx(0) = 1.d0
DO i=1,p
powerx(i) = powerx(i-1)*x
END DO
DO i=0,p-1
f(i,0) = DOT_PRODUCT(ppform(:,i),powerx(:))
END DO
f(p,0) = 1.d0 - SUM(f(0:p-1,0))
f(p,1) = - SUM(f(0:p-1,1))
END SELECT
END SUBROUTINE my_splines
!+++
END SUBROUTINE getgradr
!===========================================================================
SUBROUTINE getgradz(sp, xp, yp, f00, f10, f01)
!
! Compute the function f00 and its derivatives
! f10 = d/dx f
! f01 = d/dy f
! assuming that its PPFORM/BCOEFSC was already computed!
!
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp
DOUBLE COMPLEX, DIMENSION(:), INTENT(out) :: f00, f10, f01
!
INTEGER :: np
DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp))
INTEGER :: leftx(SIZE(xp)), lefty(SIZE(yp))
INTEGER :: i, ip, ii, jj, nidbas(2)
DOUBLE COMPLEX :: temp0(SIZE(xp),sp%sp2%order), temp1(SIZE(xp),sp%sp2%order)
DOUBLE PRECISION, ALLOCATABLE, SAVE :: funx(:,:), funy(:,:)
DOUBLE COMPLEX, ALLOCATABLE, SAVE :: ftemp0(:), ftemp1(:)
LOGICAL :: nlppform
!
! Apply periodicity if required
!
np = SIZE(xp)
nidbas(1) = sp%sp1%order-1
nidbas(2) = sp%sp2%order-1
nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform
!
! Locate the interval containing x, y
!
CALL locintv(sp%sp1, xp, leftx)
CALL locintv(sp%sp2, yp, lefty)
x(:) = xp(:) - sp%sp1%knots(leftx(:))
y(:) = yp(:) - sp%sp2%knots(lefty(:))
!
! Compute function/derivatives
!
IF(nlppform) THEN
!
! Using PPFORM
!----------
DO i=1,np
CALL my_ppval1(nidbas(1), x(i), sp%ppformz(:,leftx(i)+1,:,lefty(i)+1), &
& temp0(i,:), temp1(i,:))
END DO
!
CALL my_ppval0(nidbas(2), y, temp0, 0, f00)
CALL my_ppval0(nidbas(2), y, temp0, 1, f01)
CALL my_ppval0(nidbas(2), y, temp1, 0, f10)
ELSE
!
! Using spline expansion with sp%bcoefsc
!----------
IF(.NOT.ALLOCATED(funx)) THEN
ALLOCATE(funx(0:nidbas(1),0:1)) ! Spline and its first derivative
ALLOCATE(funy(0:nidbas(2),0:1))
ALLOCATE(ftemp0(0:nidbas(1)))
ALLOCATE(ftemp1(0:nidbas(1)))
END IF
!
DO ip=1,np
CALL my_splines(nidbas(1), x(ip), sp%sp1%val0(:,:,leftx(ip)+1), funx)
CALL my_splines(nidbas(2), y(ip), sp%sp2%val0(:,:,lefty(ip)+1), funy)
DO ii=0,nidbas(1)
ftemp0(ii) = (0.d0,0.d0)
ftemp1(ii) = (0.d0,0.d0)
DO jj=0,nidbas(2)
ftemp0(ii) = ftemp0(ii) + sp%bcoefsc(leftx(ip)+1+ii,lefty(ip)+1+jj)*funy(jj,0)
ftemp1(ii) = ftemp1(ii) + sp%bcoefsc(leftx(ip)+1+ii,lefty(ip)+1+jj)*funy(jj,1)
END DO
END DO
f00(ip) = SUM(funx(:,0)*ftemp0(:))
f01(ip) = SUM(funx(:,0)*ftemp1(:))
f10(ip) = SUM(funx(:,1)*ftemp0(:))
END DO
!-----------
END IF
CONTAINS
!+++
SUBROUTINE my_ppval0(p, x, ppform, jder, f)
!
! Compute function and derivatives from the PP representation
! for many points x(:)
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE COMPLEX, INTENT(in) :: ppform(:,:)
INTEGER, INTENT(in) :: jder
DOUBLE COMPLEX, INTENT(out) :: f(:)
DOUBLE PRECISION :: fact
INTEGER :: j
SELECT CASE (jder)
CASE(0) ! function value
SELECT CASE(p)
CASE(1)
f(:) = ppform(:,1) + x(:)*ppform(:,2)
CASE(2)
f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*ppform(:,3))
CASE(3)
f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*(ppform(:,3)+x(:)*ppform(:,4)))
CASE(4:)
DO j=p+1,1,-1
f(:) = f(:)*x(:) + ppform(:,j)
END DO
END SELECT
CASE(1) ! 1st derivative
SELECT CASE(p)
CASE(1)
f(:) = ppform(:,2)
CASE(2)
f(:) = ppform(:,2) + x(:)*2.d0*ppform(:,3)
CASE(3)
f(:) = ppform(:,2) + x(:)*(2.d0*ppform(:,3)+x(:)*3.0d0*ppform(:,4))
CASE(4:)
DO j=p,1,-1
f(:) = f(:)*x(:) + j*ppform(:,j+1)
END DO
END SELECT
CASE default ! 2nd and higher derivatives
f(:) = ppform(:,p+1)
fact = p-jder
DO j=p,jder+1,-1
f(:) = f(:)/fact*j*x(:) + ppform(:,j)
fact = fact-1.0d0
END DO
DO j=2,jder
f(:) = f(:)*j
END DO
END SELECT
END SUBROUTINE my_ppval0
!+++
SUBROUTINE my_ppval1(p, x, ppform, f0, f1)
!
! Compute function and first derivative from the PP representation
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE COMPLEX, INTENT(in) :: ppform(:,:)
DOUBLE COMPLEX, INTENT(out) :: f0(:)
DOUBLE COMPLEX, INTENT(out) :: f1(:)
DOUBLE PRECISION :: fact
INTEGER :: j
SELECT CASE(p)
CASE(1)
f0(:) = ppform(1,:) + x*ppform(2,:)
f1(:) = ppform(2,:)
CASE(2)
f0(:) = ppform(1,:) + x*(ppform(2,:)+x*ppform(3,:))
f1(:) = ppform(2,:) + x*2.d0*ppform(3,:)
CASE(3)
f0(:) = ppform(1,:) + x*(ppform(2,:)+x*(ppform(3,:)+x*ppform(4,:)))
f1(:) = ppform(2,:) + x*(2.d0*ppform(3,:)+x*3.0d0*ppform(4,:))
CASE(4:)
f0 = ppform(p+1,:)
f1 = f0
DO j=p,2,-1
f0(:) = ppform(j,:) + x*f0(:)
f1(:) = f0(:) + x*f1(:)
END DO
f0(:) = ppform(1,:) + x*f0(:)
END SELECT
END SUBROUTINE my_ppval1
!+++
SUBROUTINE my_splines(p, x, ppform, f)
INTEGER, INTENT(in) :: p
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE PRECISION, INTENT(in) :: ppform(0:p,0:p)
DOUBLE PRECISION, INTENT(out) :: f(0:p,0:1)
INTEGER :: i
DOUBLE PRECISION :: powerx(0:p)
SELECT CASE(p)
CASE(1)
f(0,0) = ppform(0,0) + x*ppform(1,0)
f(0,1) = ppform(1,0)
f(1,0) = 1.0-f(0,0)
f(1,1) = -f(0,1)
CASE(2)
f(0,0) = ppform(0,0) + x*(ppform(1,0)+x*ppform(2,0))
f(0,1) = ppform(1,0) + 2.d0*x*ppform(2,0)
f(1,0) = ppform(0,1) + x*(ppform(1,1)+x*ppform(2,1))
f(1,1) = ppform(1,1) + 2.d0*x*ppform(2,1)
f(2,0) = 1.0 - f(0,0) - f(1,0)
f(2,1) = - f(0,1) - f(1,1)
CASE(3)
f(0,0) = ppform(0,0) + x*(ppform(1,0)+x*(ppform(2,0)+x*ppform(3,0)))
f(0,1) = ppform(1,0) + x*(2.d0*ppform(2,0)+3.d0*x*ppform(3,0))
f(1,0) = ppform(0,1) + x*(ppform(1,1)+x*(ppform(2,1)+x*ppform(3,1)))
f(1,1) = ppform(1,1) + x*(2.d0*ppform(2,1)+3.d0*x*ppform(3,1))
f(2,0) = ppform(0,2) + x*(ppform(1,2)+x*(ppform(2,2)+x*ppform(3,2)))
f(2,1) = ppform(1,2) + x*(2.d0*ppform(2,2)+3.d0*x*ppform(3,2))
f(3,0) = 1.0 - f(0,0) - f(1,0) - f(2,0)
f(3,1) = - f(0,1) - f(1,1) - f(2,1)
CASE(4:)
powerx(0) = 1.d0
DO i=1,p
powerx(i) = powerx(i-1)*x
END DO
DO i=0,p-1
f(i,0) = DOT_PRODUCT(ppform(:,i),powerx(:))
END DO
f(p,0) = 1.d0 - SUM(f(0:p-1,0))
f(p,1) = - SUM(f(0:p-1,1))
END SELECT
END SUBROUTINE my_splines
!+++
END SUBROUTINE getgradz
END MODULE bsplines