Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/bsplpp.f90
c723281d03b6feature/refactoring

Last Change
Raw File

bsplpp.f90
View Options

!>
!> @file bsplpp.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 bsplpp ( t, bcoef, n, k, scrtch, break, coef, l )
!*************************************************************************
!
!! BSPLPP converts from B-spline to piecewise polynomial form.
!
! Discussion:
!
! The B-spline representation of a spline is ( T, BCOEF, N, K ),
! while the piecewise polynomial representation is
! ( BREAK, COEF, L, K ).
!
! For each breakpoint interval, the K relevant B-spline coefficients
! of the spline are found and then differenced repeatedly to get the
! B-spline coefficients of all the derivatives of the spline on that
! interval.
!
! The spline and its first K-1 derivatives are then evaluated at the
! left end point of that interval, using BSPLVB repeatedly to obtain
! the values of all B-splines of the appropriate order at that point.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) T(N+K), the knot sequence.
!
! Input, real ( kind = 8 ) BCOEF(N), the B spline coefficient sequence.
!
! Input, integer N, the number of B spline coefficients.
!
! Input, integer K, the order of the spline.
!
! Work array, real ( kind = 8 ) SCRTCH(K,K).
!
! Output, real ( kind = 8 ) BREAK(L+1), the piecewise polynomial breakpoint
! sequence. BREAK contains the distinct points in the
! sequence T(K),...,T(N+1)
!
! Output, real ( kind = 8 ) COEF(K,N), with COEF(I,J) = (I-1)st derivative
! of the spline at BREAK(J) from the right.
!
! Output, integer L, the number of polynomial pieces which
! make up the spline in the interval (T(K),T(N+1)).
!
implicit none
integer k
integer l
integer n
real ( kind = 8 ) bcoef(n)
real ( kind = 8 ) biatx(k)
real ( kind = 8 ) break(*)
real ( kind = 8 ) coef(k,n)
real ( kind = 8 ) diff
integer i
integer j
integer jp1
integer left
integer lsofar
real ( kind = 8 ) scrtch(k,k)
real ( kind = 8 ) sum1
real ( kind = 8 ) t(n+k)
lsofar = 0
break(1) = t(k)
do left = k, n
!
! Find the next nontrivial knot interval.
!
if ( t(left+1) == t(left) ) then
cycle
end if
lsofar = lsofar + 1
break(lsofar+1) = t(left+1)
if ( k <= 1 ) then
coef(1,lsofar) = bcoef(left)
cycle
end if
!
! Store the K B-spline coefficients relevant to current knot
! interval in SCRTCH(*,1).
!
do i = 1, k
scrtch(i,1) = bcoef(left-k+i)
end do
!
! For j=1,...,k-1, compute the k-j b-spline coefficients relevant to
! current knot interval for the j-th derivative by differencing
! those for the (j-1)st derivative, and store in scrtch(.,j+1) .
!
do jp1 = 2, k
j = jp1-1
do i = 1, k-j
diff = t(left+i)-t(left+i-(k-j))
if ( 0.0D+00 < diff ) then
scrtch(i,jp1)=((scrtch(i+1,j)-scrtch(i,j)) / diff ) &
* real ( k - j, kind = 8 )
end if
end do
end do
!
! For J=0, ..., K-1, find the values at T(left) of the j+1
! B-splines of order J+1 whose support contains the current
! knot interval from those of order J (in biatx ), then comb-
! ine with the B-spline coefficients (in scrtch(.,k-j) ) found earlier
! to compute the (k-j-1)st derivative at t(left) of the given
! spline.
!
call bsplvb ( t, 1, 1, t(left), left, biatx )
coef(k,lsofar) = scrtch(1,k)
do jp1 = 2, k
call bsplvb ( t, jp1, 2, t(left), left, biatx )
sum1 = 0.0D+00
do i = 1, jp1
sum1 = sum1 + biatx(i) * scrtch(i,k+1-jp1)
end do
coef(k+1-jp1,lsofar) = sum1
end do
end do
l = lsofar
return
end
c4science · Help