Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/splint.f90
c723281d03b6feature/refactoring

Last Change
Raw File

splint.f90
View Options

!>
!> @file splint.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 splint ( tau, gtau, t, n, k, q, bcoef, iflag )
!*************************************************************************
!
!! SPLINT produces the B-spline coefficients BCOEF of an interpolating spline.
!
! Discussion:
!
! The spline is of order K with knots T(1:N+K), and takes on the
! value GTAU(I) at TAU(I), for I = 1 to N.
!
! The I-th equation of the linear system
! A * BCOEF = B
! for the b-coefficients of the interpolant enforces interpolation
! at TAU(1:N).
!
! Hence, b(i)=gtau(i), all i, and a is a band matrix with 2k-1
! bands (if it is invertible).
!
! The matrix A is generated row by row and stored, diagonal by di-
! agonal, in the rows of the array q , with the main diagonal go-
! ing into row K. see comments in the program below.
!
! The banded system is then solved by a call to banfac (which con-
! structs the triangular factorization for a and stores it again in
! q ), followed by a call to banslv (which then obtains the solution
! bcoef by substitution).
!
! BANFAC does no pivoting, since the total positivity of the matrix
! A makes this unnecessary.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) TAU(N), the data point abscissas. The entries in
! TAU should be strictly increasing.
!
! Input, real ( kind = 8 ) GTAU(N), the data ordinates.
!
! Input, real ( kind = 8 ) T(N+K), the knot sequence.
!
! Input, integer N, the number of data points.
!
! Input, integer K, the order of the spline.
!
! output
!
! q, array of size (2*k-1)*n , containing the triangular factoriz-
! ation of the coefficient matrix of the linear system for the b-
! coefficients of the spline interpolant.
! the b-coefficients for the interpolant of an additional data set can
! be obtained without going through all the calculations in this
! routine, simply by loading htau into bcoef and then execut-
! ing the call banslv ( q, 2*k-1, n, k-1, k-1, bcoef )
!
! bcoef, the b-coefficients of the interpolant, of length n.
!
! iflag, an integer indicating success (= 1) or failure (= 2)
! the linear system to be solved is (theoretically) invertible if
! and only if
! t(i) < tau(i) < tau(i+k), all i.
! violation of this condition is certain to lead to iflag=2 .
!
implicit none
integer n
real ( kind = 8 ) bcoef(n)
real ( kind = 8 ) gtau(n)
integer i
integer iflag
integer ilp1mx
integer j
integer jj
integer k
integer kpkm2
integer left
real ( kind = 8 ) q((2*k-1)*n)
real ( kind = 8 ) t(n+k)
real ( kind = 8 ) tau(n)
real ( kind = 8 ) taui
kpkm2 = 2*(k-1)
left = k
do i = 1, (2*k-1)*n
q(i) = 0.0D+00
end do
!
! loop over i to construct the n interpolation equations
!
do i = 1, n
taui = tau(i)
ilp1mx = min(i+k,n+1)
!
! find left in the closed interval (i,i+k-1) such that
! t(left) <= tau(i) < t(left+1)
! matrix is singular if this is not possible
!
left = max(left,i)
if ( taui < t(left)) then
go to 70
end if
20 continue
if ( taui < t(left+1)) then
go to 30
end if
left = left+1
if ( left < ilp1mx) then
go to 20
end if
left = left-1
if ( t(left+1) < taui ) then
go to 70
end if
!
! The i-th equation enforces interpolation at taui, hence
! a(i,j)=b(j,k,t)(taui), all j. only the k entries with j =
! left-k+1,...,left actually might be nonzero. these k numbers
! are returned, in bcoef (used for temp.storage here), by the
! following
!
30 continue
call bsplvb(t,k,1,taui,left,bcoef)
!
! We therefore want bcoef(j)=b(left-k+j)(taui) to go into
! a(i,left-k+j), i.e., into q(i-(left+j)+2*k,(left+j)-k) since
! a(i+j,j) is to go into q(i+k,j), all i,j, if we consider q
! as a two-dim. array , with 2*k-1 rows (see comments in
! banfac). in the present program, we treat q as an equivalent
! one-dimensional array (because of fortran restrictions on
! dimension statements) . we therefore want bcoef(j) to go into
! entry
! i -(left+j)+2*k + ((left+j)-k-1)*(2*k-1)
! = i-left+1+(left -k)*(2*k-1) + (2*k-2)*j
! of q .
!
jj = i-left+1+(left-k)*(k+k-1)
do j = 1,k
jj = jj+kpkm2
q(jj) = bcoef(j)
end do
end do
!
! Obtain factorization of A, stored again in Q.
!
call banfac ( q, k+k-1, n, k-1, k-1, iflag )
if ( iflag == 2 ) then
write(*,*)' '
write(*,*)'SPLINT - Fatal Error!'
write(*,*)' The linear system is not invertible!'
return
end if
!
! Solve a*bcoef=gtau by backsubstitution
!
bcoef(1:n) = gtau(1:n)
call banslv ( q, k+k-1, n, k-1, k-1, bcoef )
return
70 iflag=2
write ( *, * ) ' '
write ( *, * ) 'SPLINT - Fatal Error!'
write ( *, * ) ' The linear system is not invertible!'
return
end
c4science · Help