1*DECK DEFC 2 SUBROUTINE DEFC (NDATA, XDATA, YDATA, SDDATA, NORD, NBKPT, BKPT, 3 + MDEIN, MDEOUT, COEFF, LW, W) 4C***BEGIN PROLOGUE DEFC 5C***PURPOSE Fit a piecewise polynomial curve to discrete data. 6C The piecewise polynomials are represented as B-splines. 7C The fitting is done in a weighted least squares sense. 8C***LIBRARY SLATEC 9C***CATEGORY K1A1A1, K1A2A, L8A3 10C***TYPE DOUBLE PRECISION (EFC-S, DEFC-D) 11C***KEYWORDS B-SPLINE, CONSTRAINED LEAST SQUARES, CURVE FITTING 12C***AUTHOR Hanson, R. J., (SNLA) 13C***DESCRIPTION 14C 15C This subprogram fits a piecewise polynomial curve 16C to discrete data. The piecewise polynomials are 17C represented as B-splines. 18C The fitting is done in a weighted least squares sense. 19C 20C The data can be processed in groups of modest size. 21C The size of the group is chosen by the user. This feature 22C may be necessary for purposes of using constrained curve fitting 23C with subprogram DFC( ) on a very large data set. 24C 25C For a description of the B-splines and usage instructions to 26C evaluate them, see 27C 28C C. W. de Boor, Package for Calculating with B-Splines. 29C SIAM J. Numer. Anal., p. 441, (June, 1977). 30C 31C For further discussion of (constrained) curve fitting using 32C B-splines, see 33C 34C R. J. Hanson, Constrained Least Squares Curve Fitting 35C to Discrete Data Using B-Splines, a User's 36C Guide. Sandia Labs. Tech. Rept. SAND-78-1291, 37C December, (1978). 38C 39C Input.. All TYPE REAL variables are DOUBLE PRECISION 40C NDATA,XDATA(*), 41C YDATA(*), 42C SDDATA(*) 43C The NDATA discrete (X,Y) pairs and the Y value 44C standard deviation or uncertainty, SD, are in 45C the respective arrays XDATA(*), YDATA(*), and 46C SDDATA(*). No sorting of XDATA(*) is 47C required. Any non-negative value of NDATA is 48C allowed. A negative value of NDATA is an 49C error. A zero value for any entry of 50C SDDATA(*) will weight that data point as 1. 51C Otherwise the weight of that data point is 52C the reciprocal of this entry. 53C 54C NORD,NBKPT, 55C BKPT(*) 56C The NBKPT knots of the B-spline of order NORD 57C are in the array BKPT(*). Normally the 58C problem data interval will be included between 59C the limits BKPT(NORD) and BKPT(NBKPT-NORD+1). 60C The additional end knots BKPT(I),I=1,..., 61C NORD-1 and I=NBKPT-NORD+2,...,NBKPT, are 62C required to compute the functions used to fit 63C the data. No sorting of BKPT(*) is required. 64C Internal to DEFC( ) the extreme end knots may 65C be reduced and increased respectively to 66C accommodate any data values that are exterior 67C to the given knot values. The contents of 68C BKPT(*) is not changed. 69C 70C NORD must be in the range 1 .LE. NORD .LE. 20. 71C The value of NBKPT must satisfy the condition 72C NBKPT .GE. 2*NORD. 73C Other values are considered errors. 74C 75C (The order of the spline is one more than the 76C degree of the piecewise polynomial defined on 77C each interval. This is consistent with the 78C B-spline package convention. For example, 79C NORD=4 when we are using piecewise cubics.) 80C 81C MDEIN 82C An integer flag, with one of two possible 83C values (1 or 2), that directs the subprogram 84C action with regard to new data points provided 85C by the user. 86C 87C =1 The first time that DEFC( ) has been 88C entered. There are NDATA points to process. 89C 90C =2 This is another entry to DEFC(). The sub- 91C program DEFC( ) has been entered with MDEIN=1 92C exactly once before for this problem. There 93C are NDATA new additional points to merge and 94C process with any previous points. 95C (When using DEFC( ) with MDEIN=2 it is import- 96C ant that the set of knots remain fixed at the 97C same values for all entries to DEFC( ).) 98C LW 99C The amount of working storage actually 100C allocated for the working array W(*). 101C This quantity is compared with the 102C actual amount of storage needed in DEFC( ). 103C Insufficient storage allocated for W(*) is 104C an error. This feature was included in DEFC 105C because misreading the storage formula 106C for W(*) might very well lead to subtle 107C and hard-to-find programming bugs. 108C 109C The length of the array W(*) must satisfy 110C 111C LW .GE. (NBKPT-NORD+3)*(NORD+1)+ 112C (NBKPT+1)*(NORD+1)+ 113C 2*MAX(NDATA,NBKPT)+NBKPT+NORD**2 114C 115C Output.. All TYPE REAL variables are DOUBLE PRECISION 116C MDEOUT 117C An output flag that indicates the status 118C of the curve fit. 119C 120C =-1 A usage error of DEFC( ) occurred. The 121C offending condition is noted with the SLATEC 122C library error processor, XERMSG( ). In case 123C the working array W(*) is not long enough, the 124C minimal acceptable length is printed. 125C 126C =1 The B-spline coefficients for the fitted 127C curve have been returned in array COEFF(*). 128C 129C =2 Not enough data has been processed to 130C determine the B-spline coefficients. 131C The user has one of two options. Continue 132C to process more data until a unique set 133C of coefficients is obtained, or use the 134C subprogram DFC( ) to obtain a specific 135C set of coefficients. The user should read 136C the usage instructions for DFC( ) for further 137C details if this second option is chosen. 138C COEFF(*) 139C If the output value of MDEOUT=1, this array 140C contains the unknowns obtained from the least 141C squares fitting process. These N=NBKPT-NORD 142C parameters are the B-spline coefficients. 143C For MDEOUT=2, not enough data was processed to 144C uniquely determine the B-spline coefficients. 145C In this case, and also when MDEOUT=-1, all 146C values of COEFF(*) are set to zero. 147C 148C If the user is not satisfied with the fitted 149C curve returned by DEFC( ), the constrained 150C least squares curve fitting subprogram DFC( ) 151C may be required. The work done within DEFC( ) 152C to accumulate the data can be utilized by 153C the user, if so desired. This involves 154C saving the first (NBKPT-NORD+3)*(NORD+1) 155C entries of W(*) and providing this data 156C to DFC( ) with the "old problem" designation. 157C The user should read the usage instructions 158C for subprogram DFC( ) for further details. 159C 160C Working Array.. All TYPE REAL variables are DOUBLE PRECISION 161C W(*) 162C This array is typed DOUBLE PRECISION. 163C Its length is specified as an input parameter 164C in LW as noted above. The contents of W(*) 165C must not be modified by the user between calls 166C to DEFC( ) with values of MDEIN=1,2,2,... . 167C The first (NBKPT-NORD+3)*(NORD+1) entries of 168C W(*) are acceptable as direct input to DFC( ) 169C for an "old problem" only when MDEOUT=1 or 2. 170C 171C Evaluating the 172C Fitted Curve.. 173C To evaluate derivative number IDER at XVAL, 174C use the function subprogram DBVALU( ). 175C 176C F = DBVALU(BKPT,COEFF,NBKPT-NORD,NORD,IDER, 177C XVAL,INBV,WORKB) 178C 179C The output of this subprogram will not be 180C defined unless an output value of MDEOUT=1 181C was obtained from DEFC( ), XVAL is in the data 182C interval, and IDER is nonnegative and .LT. 183C NORD. 184C 185C The first time DBVALU( ) is called, INBV=1 186C must be specified. This value of INBV is the 187C overwritten by DBVALU( ). The array WORKB(*) 188C must be of length at least 3*NORD, and must 189C not be the same as the W(*) array used in the 190C call to DEFC( ). 191C 192C DBVALU( ) expects the breakpoint array BKPT(*) 193C to be sorted. 194C 195C***REFERENCES R. J. Hanson, Constrained least squares curve fitting 196C to discrete data using B-splines, a users guide, 197C Report SAND78-1291, Sandia Laboratories, December 198C 1978. 199C***ROUTINES CALLED DEFCMN 200C***REVISION HISTORY (YYMMDD) 201C 800801 DATE WRITTEN 202C 890531 Changed all specific intrinsics to generic. (WRB) 203C 890531 REVISION DATE from Version 3.2 204C 891214 Prologue converted to Version 4.0 format. (BAB) 205C 900510 Change Prologue comments to refer to XERMSG. (RWC) 206C 900607 Editorial changes to Prologue to make Prologues for EFC, 207C DEFC, FC, and DFC look as much the same as possible. (RWC) 208C 920501 Reformatted the REFERENCES section. (WRB) 209C***END PROLOGUE DEFC 210C 211C SUBROUTINE FUNCTION/REMARKS 212C 213C DBSPVN( ) COMPUTE FUNCTION VALUES OF B-SPLINES. FROM 214C THE B-SPLINE PACKAGE OF DE BOOR NOTED ABOVE. 215C 216C DBNDAC( ), BANDED LEAST SQUARES MATRIX PROCESSORS. 217C DBNDSL( ) FROM LAWSON-HANSON, SOLVING LEAST 218C SQUARES PROBLEMS. 219C 220C DSORT( ) DATA SORTING SUBROUTINE, FROM THE 221C SANDIA MATH. LIBRARY, SAND77-1441. 222C 223C XERMSG( ) ERROR HANDLING ROUTINE 224C FOR THE SLATEC MATH. LIBRARY. 225C SEE SAND78-1189, BY R. E. JONES. 226C 227C DCOPY( ),DSCAL( ) SUBROUTINES FROM THE BLAS PACKAGE. 228C 229C WRITTEN BY R. HANSON, SANDIA NATL. LABS., 230C ALB., N. M., AUGUST-SEPTEMBER, 1980. 231C 232 DOUBLE PRECISION BKPT(*),COEFF(*),W(*),SDDATA(*),XDATA(*),YDATA(*) 233 INTEGER LW, MDEIN, MDEOUT, NBKPT, NDATA, NORD 234C 235 EXTERNAL DEFCMN 236C 237 INTEGER LBF, LBKPT, LG, LPTEMP, LWW, LXTEMP, MDG, MDW 238C 239C***FIRST EXECUTABLE STATEMENT DEFC 240C LWW=1 USAGE IN DEFCMN( ) OF W(*).. 241C LWW,...,LG-1 W(*,*) 242C 243C LG,...,LXTEMP-1 G(*,*) 244C 245C LXTEMP,...,LPTEMP-1 XTEMP(*) 246C 247C LPTEMP,...,LBKPT-1 PTEMP(*) 248C 249C LBKPT,...,LBF BKPT(*) (LOCAL TO DEFCMN( )) 250C 251C LBF,...,LBF+NORD**2 BF(*,*) 252C 253 MDG = NBKPT+1 254 MDW = NBKPT-NORD+3 255 LWW = 1 256 LG = LWW + MDW*(NORD+1) 257 LXTEMP = LG + MDG*(NORD+1) 258 LPTEMP = LXTEMP + MAX(NDATA,NBKPT) 259 LBKPT = LPTEMP + MAX(NDATA,NBKPT) 260 LBF = LBKPT + NBKPT 261 CALL DEFCMN(NDATA,XDATA,YDATA,SDDATA, 262 1 NORD,NBKPT,BKPT, 263 2 MDEIN,MDEOUT, 264 3 COEFF, 265 4 W(LBF),W(LXTEMP),W(LPTEMP),W(LBKPT), 266 5 W(LG),MDG,W(LWW),MDW,LW) 267 RETURN 268 END 269