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