relibc_openlibm/slatec/defc.f

*DECK DEFC
      SUBROUTINE DEFC (NDATA, XDATA, YDATA, SDDATA, NORD, NBKPT, BKPT,
     +   MDEIN, MDEOUT, COEFF, LW, W)
C***BEGIN PROLOGUE  DEFC
C***PURPOSE  Fit a piecewise polynomial curve to discrete data.
C            The piecewise polynomials are represented as B-splines.
C            The fitting is done in a weighted least squares sense.
C***LIBRARY   SLATEC
C***CATEGORY  K1A1A1, K1A2A, L8A3
C***TYPE      DOUBLE PRECISION (EFC-S, DEFC-D)
C***KEYWORDS  B-SPLINE, CONSTRAINED LEAST SQUARES, CURVE FITTING
C***AUTHOR  Hanson, R. J., (SNLA)
C***DESCRIPTION
C
C      This subprogram fits a piecewise polynomial curve
C      to discrete data.  The piecewise polynomials are
C      represented as B-splines.
C      The fitting is done in a weighted least squares sense.
C
C      The data can be processed in groups of modest size.
C      The size of the group is chosen by the user.  This feature
C      may be necessary for purposes of using constrained curve fitting
C      with subprogram DFC( ) on a very large data set.
C
C      For a description of the B-splines and usage instructions to
C      evaluate them, see
C
C      C. W. de Boor, Package for Calculating with B-Splines.
C                     SIAM J. Numer. Anal., p. 441, (June, 1977).
C
C      For further discussion of (constrained) curve fitting using
C      B-splines, see
C
C      R. J. Hanson, Constrained Least Squares Curve Fitting
C                   to Discrete Data Using B-Splines, a User's
C                   Guide. Sandia Labs. Tech. Rept. SAND-78-1291,
C                   December, (1978).
C
C  Input.. All TYPE REAL variables are DOUBLE PRECISION
C      NDATA,XDATA(*),
C      YDATA(*),
C      SDDATA(*)
C                         The NDATA discrete (X,Y) pairs and the Y value
C                         standard deviation or uncertainty, SD, are in
C                         the respective arrays XDATA(*), YDATA(*), and
C                         SDDATA(*).  No sorting of XDATA(*) is
C                         required.  Any non-negative value of NDATA is
C                         allowed.  A negative value of NDATA is an
C                         error.  A zero value for any entry of
C                         SDDATA(*) will weight that data point as 1.
C                         Otherwise the weight of that data point is
C                         the reciprocal of this entry.
C
C      NORD,NBKPT,
C      BKPT(*)
C                         The NBKPT knots of the B-spline of order NORD
C                         are in the array BKPT(*).  Normally the
C                         problem data interval will be included between
C                         the limits BKPT(NORD) and BKPT(NBKPT-NORD+1).
C                         The additional end knots BKPT(I),I=1,...,
C                         NORD-1 and I=NBKPT-NORD+2,...,NBKPT, are
C                         required to compute the functions used to fit
C                         the data.  No sorting of BKPT(*) is required.
C                         Internal to DEFC( ) the extreme end knots may
C                         be reduced and increased respectively to
C                         accommodate any data values that are exterior
C                         to the given knot values.  The contents of
C                         BKPT(*) is not changed.
C
C                         NORD must be in the range 1 .LE. NORD .LE. 20.
C                         The value of NBKPT must satisfy the condition
C                         NBKPT .GE. 2*NORD.
C                         Other values are considered errors.
C
C                         (The order of the spline is one more than the
C                         degree of the piecewise polynomial defined on
C                         each interval.  This is consistent with the
C                         B-spline package convention.  For example,
C                         NORD=4 when we are using piecewise cubics.)
C
C        MDEIN
C                         An integer flag, with one of two possible
C                         values (1 or 2), that directs the subprogram
C                         action with regard to new data points provided
C                         by the user.
C
C                         =1  The first time that DEFC( ) has been
C                         entered.  There are NDATA points to process.
C
C                         =2  This is another entry to DEFC().  The sub-
C                         program DEFC( ) has been entered with MDEIN=1
C                         exactly once before for this problem.  There
C                         are NDATA new additional points to merge and
C                         process with any previous points.
C                         (When using DEFC( ) with MDEIN=2 it is import-
C                         ant that the set of knots remain fixed at the
C                         same values for all entries to DEFC( ).)
C       LW
C                         The amount of working storage actually
C                         allocated for the working array W(*).
C                         This quantity is compared with the
C                         actual amount of storage needed in DEFC( ).
C                         Insufficient storage allocated for W(*) is
C                         an error.  This feature was included in DEFC
C                         because misreading the storage formula
C                         for W(*) might very well lead to subtle
C                         and hard-to-find programming bugs.
C
C                         The length of the array W(*) must satisfy
C
C                         LW .GE. (NBKPT-NORD+3)*(NORD+1)+
C                                 (NBKPT+1)*(NORD+1)+
C                               2*MAX(NDATA,NBKPT)+NBKPT+NORD**2
C
C  Output.. All TYPE REAL variables are DOUBLE PRECISION
C      MDEOUT
C                         An output flag that indicates the status
C                         of the curve fit.
C
C                         =-1  A usage error of DEFC( ) occurred.  The
C                         offending condition is noted with the SLATEC
C                         library error processor, XERMSG( ).  In case
C                         the working array W(*) is not long enough, the
C                         minimal acceptable length is printed.
C
C                         =1  The B-spline coefficients for the fitted
C                         curve have been returned in array COEFF(*).
C
C                         =2  Not enough data has been processed to
C                         determine the B-spline coefficients.
C                         The user has one of two options.  Continue
C                         to process more data until a unique set
C                         of coefficients is obtained, or use the
C                         subprogram DFC( ) to obtain a specific
C                         set of coefficients.  The user should read
C                         the usage instructions for DFC( ) for further
C                         details if this second option is chosen.
C      COEFF(*)
C                         If the output value of MDEOUT=1, this array
C                         contains the unknowns obtained from the least
C                         squares fitting process.  These N=NBKPT-NORD
C                         parameters are the B-spline coefficients.
C                         For MDEOUT=2, not enough data was processed to
C                         uniquely determine the B-spline coefficients.
C                         In this case, and also when MDEOUT=-1, all
C                         values of COEFF(*) are set to zero.
C
C                         If the user is not satisfied with the fitted
C                         curve returned by DEFC( ), the constrained
C                         least squares curve fitting subprogram DFC( )
C                         may be required.  The work done within DEFC( )
C                         to accumulate the data can be utilized by
C                         the user, if so desired.  This involves
C                         saving the first (NBKPT-NORD+3)*(NORD+1)
C                         entries of W(*) and providing this data
C                         to DFC( ) with the "old problem" designation.
C                         The user should read the usage instructions
C                         for subprogram DFC( ) for further details.
C
C  Working Array.. All TYPE REAL variables are DOUBLE PRECISION
C      W(*)
C                         This array is typed DOUBLE PRECISION.
C                         Its length is  specified as an input parameter
C                         in LW as noted above.  The contents of W(*)
C                         must not be modified by the user between calls
C                         to DEFC( ) with values of MDEIN=1,2,2,... .
C                         The first (NBKPT-NORD+3)*(NORD+1) entries of
C                         W(*) are acceptable as direct input to DFC( )
C                         for an "old problem" only when MDEOUT=1 or 2.
C
C  Evaluating the
C  Fitted Curve..
C                         To evaluate derivative number IDER at XVAL,
C                         use the function subprogram DBVALU( ).
C
C                         F = DBVALU(BKPT,COEFF,NBKPT-NORD,NORD,IDER,
C                                      XVAL,INBV,WORKB)
C
C                         The output of this subprogram will not be
C                         defined unless an output value of MDEOUT=1
C                         was obtained from DEFC( ), XVAL is in the data
C                         interval, and IDER is nonnegative and .LT.
C                         NORD.
C
C                         The first time DBVALU( ) is called, INBV=1
C                         must be specified.  This value of INBV is the
C                         overwritten by DBVALU( ).  The array WORKB(*)
C                         must be of length at least 3*NORD, and must
C                         not be the same as the W(*) array used in the
C                         call to DEFC( ).
C
C                         DBVALU( ) expects the breakpoint array BKPT(*)
C                         to be sorted.
C
C***REFERENCES  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***ROUTINES CALLED  DEFCMN
C***REVISION HISTORY  (YYMMDD)
C   800801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900510  Change Prologue comments to refer to XERMSG.  (RWC)
C   900607  Editorial changes to Prologue to make Prologues for EFC,
C           DEFC, FC, and DFC look as much the same as possible.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DEFC
C
C      SUBROUTINE           FUNCTION/REMARKS
C
C      DBSPVN( )          COMPUTE FUNCTION VALUES OF B-SPLINES.  FROM
C                         THE B-SPLINE PACKAGE OF DE BOOR NOTED ABOVE.
C
C      DBNDAC( ),         BANDED LEAST SQUARES MATRIX PROCESSORS.
C      DBNDSL( )          FROM LAWSON-HANSON, SOLVING LEAST
C                         SQUARES PROBLEMS.
C
C      DSORT( )           DATA SORTING SUBROUTINE, FROM THE
C                         SANDIA MATH. LIBRARY, SAND77-1441.
C
C      XERMSG( )          ERROR HANDLING ROUTINE
C                         FOR THE SLATEC MATH. LIBRARY.
C                         SEE SAND78-1189, BY R. E. JONES.
C
C      DCOPY( ),DSCAL( )  SUBROUTINES FROM THE BLAS PACKAGE.
C
C                         WRITTEN BY R. HANSON, SANDIA NATL. LABS.,
C                         ALB., N. M., AUGUST-SEPTEMBER, 1980.
C
      DOUBLE PRECISION BKPT(*),COEFF(*),W(*),SDDATA(*),XDATA(*),YDATA(*)
      INTEGER LW, MDEIN, MDEOUT, NBKPT, NDATA, NORD
C
      EXTERNAL DEFCMN
C
      INTEGER LBF, LBKPT, LG, LPTEMP, LWW, LXTEMP, MDG, MDW
C
C***FIRST EXECUTABLE STATEMENT  DEFC
C     LWW=1               USAGE IN DEFCMN( ) OF W(*)..
C     LWW,...,LG-1        W(*,*)
C
C     LG,...,LXTEMP-1     G(*,*)
C
C     LXTEMP,...,LPTEMP-1 XTEMP(*)
C
C     LPTEMP,...,LBKPT-1  PTEMP(*)
C
C     LBKPT,...,LBF       BKPT(*) (LOCAL TO DEFCMN( ))
C
C     LBF,...,LBF+NORD**2 BF(*,*)
C
      MDG = NBKPT+1
      MDW = NBKPT-NORD+3
      LWW = 1
      LG = LWW + MDW*(NORD+1)
      LXTEMP = LG + MDG*(NORD+1)
      LPTEMP = LXTEMP + MAX(NDATA,NBKPT)
      LBKPT  = LPTEMP + MAX(NDATA,NBKPT)
      LBF = LBKPT + NBKPT
      CALL DEFCMN(NDATA,XDATA,YDATA,SDDATA,
     1         NORD,NBKPT,BKPT,
     2         MDEIN,MDEOUT,
     3         COEFF,
     4         W(LBF),W(LXTEMP),W(LPTEMP),W(LBKPT),
     5         W(LG),MDG,W(LWW),MDW,LW)
      RETURN
      END