relibc_openlibm/slatec/bspdoc.f

*DECK BSPDOC
      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
C***FIRST EXECUTABLE STATEMENT  BSPDOC
      RETURN
      END