Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/pppack/spli2d.f90
93a1e9712d04feature/refactoring

Last Change
Raw File

spli2d.f90
View Options

!>
!> @file spli2d.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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
subroutine spli2d ( tau, gtau, t, n, k, m, work, q, bcoef, iflag )
!*************************************************************************
!
!! SPLI2D produces a interpolatory tensor product spline.
!
! Discussion:
!
! SPLI2D is an extended version of SPLINT.
!
! SPLI2D produces the B-spline coefficients BCOEF(J,.) of the
! spline of order K with knots T(I), I=1,..., N+K, which takes on
! the value GTAU(I,J) at TAU(I), I=1,..., N, J=1,...,M.
!
! The I-th equation of the linear system
!
! A * BCOEF = B
!
! for the B-spline coefficients of the interpolant enforces
! interpolation at TAU(I), I=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
! diagonal, in the rows of the array Q, with the main diagonal
! going into row K.
!
! The banded system is then solved by a call to BANFAC, which
! constructs 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.
!
! Reference:
!
! Carl DeBoor,
! A Practical Guide to Splines,
! Springer Verlag.
!
! Parameters:
!
! Input, real ( kind = 8 ) TAU(N), contains the data point abscissas.
! TAU must be strictly increasing
!
! Input, real ( kind = 8 ) GTAU(N), contains the data point ordinates,
! J=1,...,M.
!
! Input, real ( kind = 8 ) T(N+K), the knot sequence.
!
! Input, integer N, the number of data points and the
! dimension of the spline space SPLINE(K,T)
!
! Input, integer K, the order of the spline.
!
! Input, integer M, the number of data sets.
!
! Work space, real ( kind = 8 ) WORK(N).
!
! Output, real ( kind = 8 ) Q(2*K-1)*N, containing the triangular
! factorization of the coefficient matrix of the linear
! system for the B-spline coefficients of the spline interpolant.
!
! The B-spline coefficients for the interpolant of an additional
! data set (TAU(I),HTAU(I)), I=1,...,N with the same data
! abscissae can be obtained without going through all the
! calculations in this routine, simply by loading HTAU into
! BCOEF and then using the statement
!
! CALL BANSLV(Q,2*K-1,N,K-1,K-1,BCOEF)
!
! Output, real ( kind = 8 ) BCOEF(N), the B-spline coefficients of
! the interpolant.
!
! Output, integer IFLAG, error indicator.
! 1, no error.
! 2, an error occurred, which may have been caused by
! singularity of the linear system.
!
! The linear system to be solved is theoretically invertible if
! and only if
!
! T(I) < TAU(I) < TAU(I+K), for all I.
!
! Violation of this condition is certain to lead to IFLAG=2.
!
implicit none
integer m
integer n
real ( kind = 8 ) bcoef(m,n)
real ( kind = 8 ) gtau(n,m)
integer i
integer iflag
integer ilp1mx
integer j
integer jj
integer k
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
real ( kind = 8 ) work(n)
left = k
do i = 1, (2*k-1)*n
q(i) = 0.0
end do
!
! Construct the N interpolation equations.
!
do i = 1, n
taui = tau(i)
ilp1mx = min(i+k,n+1)
!
! Find the index LEFT in the closed interval (I,I+K-1) such
! that:
!
! T(LEFT) < = TAU(I) < T(LEFT+1)
!
! The matrix will be singular if this is not possible.
!
left = max(left,i)
if ( taui < t(left) ) then
iflag = 2
write(*,*)' '
write(*,*)'SPLI2D - Fatal error!'
write(*,*)' The TAU array is not strictly increasing.'
stop
end if
20 continue
if ( t(left+1) <= taui ) then
left = left+1
if ( left < ilp1mx ) then
go to 20
end if
left = left-1
if ( t(left+1) < taui ) then
iflag = 2
write(*,*)' '
write(*,*)'SPLI2D - Fatal error!'
write(*,*)' The TAU array is not strictly increasing.'
stop
end if
end if
!
! The I-th equation enforces interpolation at TAUI, hence
!
! A(I,J)=B(J,K,T)(TAUI), for all J.
!
! Only the K entries with J=LEFT-K+1, ..., LEFT actually might be
! nonzero. These K numbers are returned, in WORK (used for
! temporary storage here), by the following call:
!
call bsplvb(t,k,1,taui,left,work)
!
! We therefore want
!
! WORK(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
! ?? LOST LINE ??
! 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+k+k-2
q(jj) = work(j)
end do
end do
!
! Factor A, stored again in Q.
!
call banfac(q,k+k-1,n,k-1,k-1,iflag)
if ( iflag == 2 ) then
write(*,*)' '
write(*,*)'SPLI2D - Fatal error!'
write(*,*)' BANFAC reports that the matrix is singular.'
stop
end if
!
! Solve A*BCOEF=GTAU by backsubstitution.
!
do j = 1, m
work(1:n) = gtau(1:n,j)
call banslv ( q, k+k-1, n, k-1, k-1, work )
bcoef(j,1:n) = work(1:n)
end do
return
end
c4science · Help