Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/bsplvb.f90
8f478ec0accefeature/refactoring

Last Change
Raw File

bsplvb.f90
View Options

!>
!> @file bsplvb.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 bsplvb ( t, jhigh, index, x, left, biatx )
!*************************************************************************
!
!! BSPLVB evaluates B-splines at a point X with a given knot sequence.
!
! Discusion:
!
! BSPLVB evaluates all possibly nonzero B-splines at X of order
!
! JOUT = MAX ( JHIGH, (J+1)*(INDEX-1) )
!
! with knot sequence T.
!
! The recurrence relation
!
! X - T(I) T(I+J+1) - X
! B(I,J+1)(X) = ----------- * B(I,J)(X) + --------------- * B(I+1,J)(X)
! T(I+J)-T(I) T(I+J+1)-T(I+1)
!
! is used to generate B(LEFT-J:LEFT,J+1)(X) from B(LEFT-J+1:LEFT,J)(X)
! storing the new values in BIATX over the old.
!
! The facts that
!
! B(I,1)(X) = 1 if T(I) <= X < T(I+1)
!
! and that
!
! B(I,J)(X) = 0 unless T(I) <= X < T(I+J)
!
! are used.
!
! The particular organization of the calculations follows
! algorithm 8 in chapter X of the text.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) T(LEFT+JOUT), the knot sequence. T is assumed to
! be nondecreasing, and also, T(LEFT) must be strictly less than
! T(LEFT+1).
!
! Input, integer JHIGH, INDEX, determine the order
! JOUT = MAX(JHIGH,(J+1)*(INDEX-1))
! of the B-splines whose values at X are to be returned.
! INDEX is used to avoid recalculations when several
! columns of the triangular array of B-spline values are
! needed, for example, in BVALUE or in BSPLVD.
!
! If INDEX = 1, the calculation starts from scratch and the entire
! triangular array of B-spline values of orders
! 1, 2, ...,JHIGH is generated order by order, i.e.,
! column by column.
!
! If INDEX = 2, only the B-spline values of order J+1, J+2, ..., JOUT
! are generated, the assumption being that BIATX, J,
! DELTAL, DELTAR are, on entry, as they were on exit
! at the previous call. In particular, if JHIGH = 0,
! then JOUT = J+1, i.e., just the next column of B-spline
! values is generated.
!
! WARNING: the restriction JOUT <= JMAX (= 20) is
! imposed arbitrarily by the dimension statement for DELTAL
! and DELTAR, but is nowhere checked for.
!
! Input, real ( kind = 8 ) X, the point at which the B-splines
! are to be evaluated.
!
! Input, integer LEFT, an integer chosen so that
! T(LEFT) <= X <= T(LEFT+1).
!
! Output, real ( kind = 8 ) BIATX(JOUT), with BIATX(I) containing the
! value at X of the polynomial of order JOUT which agrees
! with the B-spline B(LEFT-JOUT+I,JOUT,T) on the interval
! (T(LEFT),T(LEFT+1)).
!
implicit none
integer, parameter :: jmax = 20
integer jhigh
real ( kind = 8 ) biatx(jhigh)
!!$ real ( kind = 8 ), save, dimension ( jmax ) :: deltal
!!$ real ( kind = 8 ), save, dimension ( jmax ) :: deltar
real ( kind = 8 ), dimension ( jmax ) :: deltal
real ( kind = 8 ), dimension ( jmax ) :: deltar
integer i
integer index
!!$ integer, save :: j = 1
integer :: j
integer left
real ( kind = 8 ) saved
real ( kind = 8 ) t(left+jhigh)
real ( kind = 8 ) term
real ( kind = 8 ) x
! Forces starting always from scratch!
!!$ if ( index == 1 ) then
j = 1
biatx(1) = 1.0D+00
if ( jhigh <= j ) then
return
end if
!!$ end if
if ( t(left+1) <= t(left) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BSPLVB - Fatal error!'
write ( *, '(a)' ) ' It is required that T(LEFT) < T(LEFT+1).'
write ( *, '(a,i6)' ) ' But LEFT = ', left
write ( *, '(a,g14.6)' ) ' T(LEFT) = ', t(left)
write ( *, '(a,g14.6)' ) ' T(LEFT+1) = ', t(left+1)
stop
end if
do
deltar(j) = t(left+j) - x
deltal(j) = x - t(left+1-j)
saved = 0.0D+00
do i = 1, j
term = biatx(i) / ( deltar(i) + deltal(j+1-i) )
biatx(i) = saved + deltar(i) * term
saved = deltal(j+1-i) * term
end do
biatx(j+1) = saved
j = j + 1
if ( jhigh <= j ) then
exit
end if
end do
return
end
c4science · Help