SUBROUTINE PDA_DBINTK (X, Y, T, N, K, BCOEF, Q, WORK, STATUS)
***BEGIN PROLOGUE PDA_DBINTK
***PURPOSE Compute the B-representation of a spline which interpolates
given data.
***LIBRARY SLATEC
***CATEGORY E1A
***TYPE DOUBLE PRECISION (BINTK-S, PDA_DBINTK-D)
***KEYWORDS B-SPLINE, DATA FITTING, INTERPOLATION
***AUTHOR Amos, D. E., (SNLA)
***DESCRIPTION
Written by Carl de Boor and modified by D. E. Amos
Abstract **** a double precision routine ****
PDA_DBINTK is the PDA_SPLINT routine of the reference.
PDA_DBINTK produces the B-spline coefficients, BCOEF, of the
B-spline of order K with knots T(I), I=1,...,N+K, which
takes on the value Y(I) at X(I), I=1,...,N. The spline or
any of its derivatives can be evaluated by calls to PDA_DBVALU.
The I-th equation of the linear system A*BCOEF = B for the
coefficients of the interpolant enforces interpolation at
X(I), I=1,...,N. Hence, B(I) = Y(I), for all I, and A is
a band matrix with 2K-1 bands if A is invertible. The matrix
A is generated row by row and stored, diagonal by diagonal,
in the rows of Q, with the main diagonal going into row K.
The banded system is then solved by a call to PDA_DBNFAC (which
constructs the triangular factorization for A and stores it
again in Q), followed by a call to PDA_DBNSLV (which then
obtains the solution BCOEF by substitution). PDA_DBNFAC does no
pivoting, since the total positivity of the matrix A makes
this unnecessary. The linear system to be solved is
(theoretically) invertible if and only if
T(I) .LT. X(I) .LT. T(I+K), for all I.
Equality is permitted on the left for I=1 and on the right
for I=N when K knots are used at X(1) or X(N). Otherwise,
violation of this condition is certain to lead to an error.
Description of Arguments
Input X,Y,T are double precision
X - vector of length N containing data point abscissa
in strictly increasing order.
Y - corresponding vector of length N containing data
point ordinates.
T - knot vector of length N+K
Since T(1),..,T(K) .LE. X(1) and T(N+1),..,T(N+K)
.GE. X(N), this leaves only N-K knots (not nec-
essarily X(I) values) interior to (X(1),X(N))
N - number of data points, N .GE. K
K - order of the spline, K .GE. 1
Output BCOEF,Q,WORK are double precision
BCOEF - a vector of length N containing the B-spline
coefficients
Q - a work vector of length (2*K-1)*N, containing
the triangular factorization of the coefficient
matrix of the linear system being solved. The
coefficients for the interpolant of an
additional data set (X(I),YY(I)), I=1,...,N
with the same abscissa can be obtained by loading
YY into BCOEF and then executing
CALL PDA_DBNSLV (Q,2K-1,N,K-1,K-1,BCOEF)
WORK - work vector of length 2*K
STATUS - Returned error status.
The status must be zero on entry. This
routine does not check the status on entry.
Error Conditions
Improper input is a fatal error
Singular system of equations is a fatal error
***REFERENCES D. E. Amos, Computation with splines and B-splines,
Report SAND78-1968, Sandia Laboratories, March 1979.
Carl de Boor, Package for calculating with B-splines,
SIAM Journal on Numerical Analysis 14, 3 (June 1977),
pp. 441-472.
Carl de Boor, A Practical Guide to Splines, Applied
Mathematics Series 27, Springer-Verlag, New York,
1978.
***ROUTINES CALLED PDA_DBNFAC, PDA_DBNSLV, PDA_DBSPVN, PDA_XERMSG
***REVISION HISTORY (YYMMDD)
800901 DATE WRITTEN
890531 Changed all specific intrinsics to generic. (WRB)
890831 Modified array declarations. (WRB)
890831 REVISION DATE from Version 3.2
891214 Prologue converted to Version 4.0 format. (BAB)
900315 CALLs to XERROR changed to CALLs to PDA_XERMSG. (THJ)
900326 Removed duplicate information from DESCRIPTION section.
(WRB)
920501 Reformatted the REFERENCES section. (WRB)
950403 Implement status. (HME)
***END PROLOGUE PDA_DBINTK
PDA [1ex