# bspdoc.f

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