*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