function bspdoc(varargin)
%***BEGIN PROLOGUE BSPDOC
%***PURPOSE Documentation for BSPLINE, a package of subprograms for
% working with piecewise polynomial functions
% in B-representation.
%***LIBRARY SLATEC
%***CATEGORY E, E1A, K, Z
%***TYPE ALL (BSPDOC-A)
%***KEYWORDS B-SPLINE, DOCUMENTATION, SPLINES
%***AUTHOR Amos, D. E., (SNLA)
%***DESCRIPTION
%
% Abstract
% BSPDOC is a non-executable, B-spline documentary routine.
% The narrative describes a B-spline and the routines
% necessary to manipulate B-splines at a fairly high level.
% The basic package described herein is that of reference
% 5 with names altered to prevent duplication and conflicts
% with routines from reference 3. The call lists used here
% are also different. Work vectors were added to ensure
% portability and proper execution in an overlay environ-
% ment. These work arrays can be used for other purposes
% except as noted in BSPVN. While most of the original
% routines in reference 5 were restricted to orders 20
% or less, this restriction was removed from all routines
% except the quadrature routine BSQAD. (See the section
% below on differentiation and integration for details.)
%
% The subroutines referenced below are single precision
% routines. Corresponding doubleprecision versions are also
% part of the package, and these are referenced by prefixing
% a D in front of the single precision name. For example,
% BVALU and DBVALU are the single and doubleprecision
% versions for evaluating a B-spline or any of its deriva-
% tives in the B-representation.
%
% ****Description of B-Splines****
%
% A collection of polynomials of fixed degree K-1 defined on a
% subdivision (X(I),X(I+1)), I=1,...,M-1 of (A,B) with X(1)=A,
% X(M)=B is called a B-spline of order K. If the spline has K-2
% continuous derivatives on (A,B), then the B-spline is simply
% called a spline of order K. Each of the M-1 polynomial pieces
% has K coefficients, making a total of K(M-1) parameters. This
% B-spline and its derivatives have M-2 jumps at the subdivision
% points X(I), I=2,...,M-1. Continuity requirements at these
% subdivision points add constraints and reduce the number of free
% parameters. If a B-spline is continuous at each of the M-2 sub-
% division points, there are K(M-1)-(M-2) free parameters; if in
% addition the B-spline has continuous first derivatives, there
% are K(M-1)-2(M-2) free parameters, etc., until we get to a
% spline where we have K(M-1)-(K-1,M-2) = M+K-2 free parameters.
% Thus, the principle is that increasing the continuity of
% derivatives decreases the number of free parameters and
% conversely.
%
% The points at which the polynomials are tied together by the
% continuity conditions are called knots. If two knots are
% allowed to come together at some X(I), then we say that we
% have a knot of multiplicity 2 there, and the knot values are
% the X(I) value. If we reverse the procedure of the first
% paragraph, we find that adding a knot to increase multiplicity
% increases the number of free parameters and, according to the
% principle above, we thereby introduce a discontinuity in what
% was the highest continuous derivative at that knot. Thus, the
% number of free parameters is N = NU+K-2 where NU is the sum
% of multiplicities at the X(I) values with X(1) and X(M) of
% multiplicity 1 (NU = M if all knots are simple, i.e., for a
% spline, all knots have multiplicity 1.) Each knot can have a
% multiplicity of at most K. A B-spline is commonly written in the
% B-representation
%
% Y(X) = sum( A(I)*B(I,X), I=1 , N)
%
% to show the explicit dependence of the spline on the free
% parameters or coefficients A(I)=BCOEF(I) and basis functions
% B(I,X). These basis functions are themselves special B-splines
% which are zero except on (at most) K adjoining intervals where
% each B(I,X) is positive and, in most cases, hat or bell-
% shaped. In order for the nonzero part of B(1,X) to be a spline
% covering (X(1),X(2)), it is necessary to put K-1 knots to the
% left of A and similarly for B(N,X) to the right of B. Thus, the
% total number of knots for this representation is NU+2K-2 = N+K.
% These knots are carried in an array T(*) dimensioned by at least
% N+K. From the construction, A=T(K) and B=T(N+1) and the spline is
% defined on T(K).LE.X.LE.T(N+1). The nonzero part of each basis
% function lies in the Interval (T(I),T(I+K)). In many problems
% where extrapolation beyond A or B is not anticipated, it is common
% practice to set T(1)=T(2)=...=T(K)=A and T(N+1)=T(N+2)=...=
% T(N+K)=B. In summary, since T(K) and T(N+1) as well as
% interior knots can have multiplicity K, the number of free
% parameters N = sum of multiplicities - K. The fact that each
% B(I,X) function is nonzero over at most K intervals means that
% for a given X value, there are at most K nonzero terms of the
% sum. This leads to banded matrices in linear algebra problems,
% and references 3 and 6 take advantage of this in con-
% structing higher level routines to achieve speed and avoid
% ill-conditioning.
%
% ****Basic Routines****
%
% The basic routines which most casual users will need are those
% concerned with direct evaluation of splines or B-splines.
% Since the B-representation, denoted by (T,BCOEF,N,K), is
% preferred because of numerical stability, the knots T(*), the
% B-spline coefficients BCOEF(*), the number of coefficients N,
% and the order K of the polynomial pieces (of degree K-1) are
% usually given. While the knot array runs from T(1) to T(N+K),
% the B-spline is normally defined on the interval T(K).LE.X.LE.
% T(N+1). To evaluate the B-spline or any of its derivatives
% on this interval, one can use
%
% Y = BVALU(T,BCOEF,N,K,ID,X,INBV,WORK)
%
% where ID is an integer for the ID-th derivative, 0.LE.ID.LE.K-1.
% ID=0 gives the zero-th derivative or B-spline value at X.
% If X.LT.T(K) or X.GT.T(N+1), whether by mistake or the result
% of round off accumulation in incrementing X, BVALU gives a
% diagnostic. INBV is an initialization parameter which is set
% to 1 on the first call. Distinct splines require distinct
% INBV parameters. WORK is a scratch vector of length at least
% 3*K.
%
% When more conventional communication is needed for publication,
% physical interpretation, etc., the B-spline coefficients can
% be converted to piecewise polynomial (PP) coefficients. Thus,
% the breakpoints (distinct knots) XI(*), the number of
% polynomial pieces LXI, and the (right) derivatives C(*,J) at
% each breakpoint XI(J) are needed to define the Taylor
% expansion to the right of XI(J) on each interval XI(J).LE.
% X.LT.XI(J+1), J=1,LXI where XI(1)=A and XI(LXI+1)=B.
% These are obtained from the (T,BCOEF,N,K) representation by
%
% CALL BSPPP(T,BCOEF,N,K,LDC,C,XI,LXI,WORK)
%
% where LDC.GE.K is the leading dimension of the matrix C and
% WORK is a scratch vector of length at least K*(N+3).
% Then the PP-representation (C,XI,LXI,K) of Y(X), denoted
% by Y(J,X) on each interval XI(J).LE.X.LT.XI(J+1), is
%
% Y(J,X) = sum( C(I,J)*((X-XI(J))**(I-1))/factorial(I-1), I=1,K)
%
% for J=1,...,LXI. One must view this conversion from the B-
% to the PP-representation with some skepticism because the
% conversion may lose significant digits when the B-spline
% varies in an almost discontinuous fashion. To evaluate
% the B-spline or any of its derivatives using the PP-
% representation, one uses
%
% Y = PPVAL(LDC,C,XI,LXI,K,ID,X,INPPV)
%
% where ID and INPPV have the same meaning and usage as ID and
% INBV in BVALU.
%
% To determine to what extent the conversion process loses
% digits, compute the relative error ABS((Y1-Y2)/Y2) over
% the X interval with Y1 from PPVAL and Y2 from BVALU. A
% major reason for considering PPVAL is that evaluation is
% much faster than that from BVALU.
%
% Recall that when multiple knots are encountered, jump type
% discontinuities in the B-spline or its derivatives occur
% at these knots, and we need to know that BVALU and PPVAL
% return right limiting values at these knots except at
% X=B where left limiting values are returned. These values
% are used for the Taylor expansions about left end points of
% breakpoint intervals. That is, the derivatives C(*,J) are
% right derivatives. Note also that a computed X value which,
% mathematically, would be a knot value may differ from the knot
% by a round off error. When this happens in evaluating a dis-
% continuous B-spline or some discontinuous derivative, the
% value at the knot and the value at X can be radically
% different. In this case, setting X to a T or XI value makes
% the computation precise. For left limiting values at knots
% other than X=B, see the prologues to BVALU and other
% routines.
%
% ****Interpolation****
%
% BINTK is used to generate B-spline parameters (T,BCOEF,N,K)
% which will interpolate the data by calls to BVALU. A similar
% interpolation can also be done for cubic splines using BINT4
% or the code in reference 7. If the PP-representation is given,
% one can evaluate this representation at an appropriate number of
% abscissas to create data then use BINTK or BINT4 to generate
% the B-representation.
%
% ****Differentiation and Integration****
%
% Derivatives of B-splines are obtained from BVALU or PPVAL.
% Integrals are obtained from BSQAD using the B-representation
% (T,BCOEF,N,K) and PPQAD using the PP-representation (C,XI,LXI,
% K). More complicated integrals involving the product of a
% of a function F and some derivative of a B-spline can be
% evaluated with BFQAD or PFQAD using the B- or PP- represen-
% tations respectively. All quadrature routines, except for PPQAD,
% are limited in accuracy to 18 digits or working precision,
% whichever is smaller. PPQAD is limited to working precision
% only. In addition, the order K for BSQAD is limited to 20 or
% less. If orders greater than 20 are required, use BFQAD with
% F(X) = 1.
%
% ****Extrapolation****
%
% Extrapolation outside the interval (A,B) can be accomplished
% easily by the PP-representation using PPVAL. However,
% caution should be exercised, especially when several knots
% are located at A or B or when the extrapolation is carried
% significantly beyond A or B. On the other hand, direct
% evaluation with BVALU outside A=T(K).LE.X.LE.T(N+1)=B
% produces an error message, and some manipulation of the knots
% and coefficients are needed to extrapolate with BVALU. This
% process is described in reference 6.
%
% ****Curve Fitting and Smoothing****
%
% Unless one has many accurate data points, direct inter-
% polation is not recommended for summarizing data. The
% results are often not in accordance with intuition since the
% fitted curve tends to oscillate through the set of points.
% Monotone splines (reference 7) can help curb this undulating
% tendency but constrained least squares is more likely to give an
% acceptable fit with fewer parameters. Subroutine FC, des-
% cribed in reference 6, is recommended for this purpose. The
% output from this fitting process is the B-representation.
%
% **** Routines in the B-Spline Package ****
%
% Single Precision Routines
%
% The subroutines referenced below are SINGLE PRECISION
% routines. Corresponding doubleprecision versions are also
% part of the package and these are referenced by prefixing
% a D in front of the single precision name. For example,
% BVALU and DBVALU are the SINGLE and doubleprecision
% versions for evaluating a B-spline or any of its deriva-
% tives in the B-representation.
%
% BINT4 - interpolates with splines of order 4
% BINTK - interpolates with splines of order k
% BSQAD - integrates the B-representation on subintervals
% PPQAD - integrates the PP-representation
% BFQAD - integrates the product of a function F and any spline
% derivative in the B-representation
% PFQAD - integrates the product of a function F and any spline
% derivative in the PP-representation
% BVALU - evaluates the B-representation or a derivative
% PPVAL - evaluates the PP-representation or a derivative
% INTRV - gets the largest index of the knot to the left of x
% BSPPP - converts from B- to PP-representation
% BSPVD - computes nonzero basis functions and derivatives at x
% BSPDR - sets up difference array for BSPEV
% BSPEV - evaluates the B-representation and derivatives
% BSPVN - called by BSPEV, BSPVD, BSPPP and BINTK for function and
% derivative evaluations
% Auxiliary Routines
%
% BSGQ8,PPGQ8,BNSLV,BNFAC,XERMSG,DBSGQ8,DPPGQ8,DBNSLV,DBNFAC
%
% Machine Dependent Routines
%
% I1MACH, R1MACH, D1MACH
%
%***REFERENCES 1. D. E. Amos, Computation with splines and
% B-splines, Report SAND78-1968, Sandia
% Laboratories, March 1979.
% 2. D. E. Amos, Quadrature subroutines for splines and
% B-splines, Report SAND79-1825, Sandia Laboratories,
% December 1979.
% 3. Carl de Boor, A Practical Guide to Splines, Applied
% Mathematics Series 27, Springer-Verlag, New York,
% 1978.
% 4. Carl de Boor, On calculating with B-Splines, Journal
% of Approximation Theory 6, (1972), pp. 50-62.
% 5. Carl de Boor, Package for calculating with B-splines,
% SIAM Journal on Numerical Analysis 14, 3 (June 1977),
% pp. 441-472.
% 6. R. J. Hanson, Constrained least squares curve fitting
% to discrete data using B-splines, a users guide,
% Report SAND78-1291, Sandia Laboratories, December
% 1978.
% 7. F. N. Fritsch and R. E. Carlson, Monotone piecewise
% cubic interpolation, SIAM Journal on Numerical Ana-
% lysis 17, 2 (April 1980), pp. 238-246.
%***ROUTINES CALLED (NONE)
%***REVISION HISTORY (YYMMDD)
% 810223 DATE WRITTEN
% 861211 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900723 PURPOSE section revised. (WRB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE BSPDOC
%***FIRST EXECUTABLE STATEMENT BSPDOC
end
%DECK BSPDR
