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