# efc.f

          SUBROUTINE EFC (NDATA, XDATA, YDATA, SDDATA, NORD, NBKPT, BKPT,
         +   MDEIN, MDEOUT, COEFF, LW, W)
    C***BEGIN PROLOGUE  EFC
    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      SINGLE PRECISION (EFC-S, DEFC-D)
    C***KEYWORDS  B-SPLINE, CURVE FITTING, WEIGHTED LEAST SQUARES
    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 FC( ) 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..
    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  EFC( ) 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 EFC( ) has been
    C                         entered.  There are NDATA points to process.
    C
    C                         =2  This is another entry to EFC( ).  The sub-
    C                         program EFC( ) 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 EFC( ) with MDEIN=2 it is import-
    C                         ant that the set of knots remain fixed at the
    C                         same values for all entries to EFC( ).)
    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 EFC( ).
    C                         Insufficient storage allocated for W(*) is
    C                         an error.  This feature was included in EFC( )
    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..
    C      MDEOUT
    C                         An output flag that indicates the status
    C                         of the curve fit.
    C
    C                         =-1  A usage error of EFC( ) 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 FC( ) to obtain a specific
    C                         set of coefficients.  The user should read
    C                         the usage instructions for FC( ) 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 EFC( ), the constrained
    C                         least squares curve fitting subprogram FC( )
    C                         may be required.  The work done within EFC( )
    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 FC( ) with the "old problem" designation.
    C                         The user should read the usage instructions
    C                         for subprogram FC( ) for further details.
    C
    C  Working Array..
    C      W(*)
    C                         This array is typed REAL.
    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 EFC( ) 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 FC( )
    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 BVALU( ).
    C
    C                         F = BVALU(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 EFC( ), XVAL is in the data
    C                         interval, and IDER is nonnegative and .LT.
    C                         NORD.
    C
    C                         The first time BVALU( ) is called, INBV=1
    C                         must be specified.  This value of INBV is the
    C                         overwritten by BVALU( ).  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 EFC( ).
    C
    C                         BVALU( ) 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  EFCMN
    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  EFC