Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/bsplvd.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

bsplvd.f90
View Options

!>
!> @file bsplvd.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine bsplvd ( t, k, x, left, a, dbiatx, nderiv )
!*************************************************************************
!
!! BSPLVD calculates the nonvanishing B-splines and derivatives at X.
!
! Discussion:
!
! Values at X of all the relevant B-splines of order K, K-1,..., K+1-NDERIV
! are generated via BSPLVB and stored temporarily in DBIATX.
!
! Then, the B-spline coefficients of the required derivatives
! of the B-splines of interest are generated by differencing,
! each from the preceding one of lower order, and combined with
! the values of B-splines of corresponding order in DBIATX
! to produce the desired values.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) T(LEFT+K), the knot sequence. It is assumed that
! T(LEFT) < T(LEFT+1). Also, the output is correct only if
! T(LEFT) <= X <= T(LEFT+1) .
!
! Input, integer K, the order of the B-splines to be evaluated.
!
! Input, real ( kind = 8 ) X, the point at which these values are sought.
!
! Input, integer LEFT, indicates the left endpoint of the interval of
! interest. The K B-splines whose support contains the interval
! (T(LEFT), T(LEFT+1)) are to be considered.
!
! Workspace, real ( kind = 8 ) A(K,K).
!
! Output, real ( kind = 8 ) DBIATX(K,NDERIV). DBIATX(I,M) contains
! the value of the (M-1)st derivative of the (LEFT-K+I)-th B-spline
! of order K for knot sequence T, I=M,...,K, M=1,...,NDERIV.
!
! Input, integer NDERIV, indicates that values of
! B-splines and their derivatives up to but not
! including the NDERIV-th are asked for.
!
implicit none
integer k
integer left
integer nderiv
real ( kind = 8 ) a(k,k)
real ( kind = 8 ) dbiatx(k,nderiv)
real ( kind = 8 ) factor
real ( kind = 8 ) fkp1mm
integer i
integer ideriv
integer il
integer j
integer jlow
integer jp1mid
integer ldummy
integer m
integer mhigh
real ( kind = 8 ) sum1
real ( kind = 8 ) t(left+k)
real ( kind = 8 ) x
mhigh = max ( min ( nderiv, k ), 1 )
!
! MHIGH is usually equal to nderiv.
!
call bsplvb ( t, k+1-mhigh, 1, x, left, dbiatx )
if ( mhigh == 1 ) then
return
end if
!
! The first column of DBIATX always contains the B-spline values
! for the current order. These are stored in column K+1-current
! order before BSPLVB is called to put values for the next
! higher order on top of it.
!
ideriv = mhigh
do m = 2, mhigh
jp1mid = 1
do j = ideriv, k
dbiatx(j,ideriv) = dbiatx(jp1mid,1)
jp1mid = jp1mid+1
end do
ideriv = ideriv-1
call bsplvb(t,k+1-ideriv,2,x,left,dbiatx)
end do
!
! At this point, b(left-k+i, k+1-j)(x) is in dbiatx(i,j) for
! i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the
! first column of dbiatx is already in final form. to obtain cor-
! ??? LOST A LINE ???
! rate their b-repr. by differencing, then evaluate at x.
!
jlow = 1
do i = 1, k
do j = jlow,k
a(j,i) = 0.0D+00
end do
jlow = i
a(i,i) = 1.0D+00
end do
!
! At this point, a(.,j) contains the b-coefficients for the J-th of the
! k b-splines of interest here.
!
do m = 2, mhigh
fkp1mm = real ( k + 1 - m, kind = 8 )
il = left
i = k
!
! For j=1,...,k, construct b-coefficients of (m-1)st derivative of
! b-splines from those for preceding derivative by differencing
! and store again in a(.,j) . The fact that a(i,j)=0 for
! i < j is used.
!
do ldummy = 1, k+1-m
factor = fkp1mm/(t(il+k+1-m)-t(il))
!
! The assumption that t(left) < t(left+1) makes denominator
! in factor nonzero.
!
do j = 1, i
a(i,j) = (a(i,j)-a(i-1,j))*factor
end do
il = il-1
i = i-1
end do
!
! For i=1,...,k, combine b-coefficients a(.,i) with B-spline values
! stored in dbiatx(.,m) to get value of (m-1)st derivative of
! i-th b-spline (of interest here) at x , and store in
! dbiatx(i,m). storage of this value over the value of a b-spline
! of order m there is safe since the remaining b-spline derivat-
! ives of the same order do not use this value due to the fact
! that a(j,i)=0 for j < i.
!
do i = 1, k
sum1 = 0.0D+00
jlow = max(i,m)
do j = jlow,k
sum1 = sum1 + a(j,i) * dbiatx(j,m)
end do
dbiatx(i,m) = sum1
end do
end do
return
end
c4science · Help