*DECK DFC SUBROUTINE DFC (NDATA, XDATA, YDATA, SDDATA, NORD, NBKPT, BKPT, + NCONST, XCONST, YCONST, NDERIV, MODE, COEFF, W, IW) C***BEGIN PROLOGUE DFC 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 Equality and inequality constraints can be imposed on the C fitted curve. C***LIBRARY SLATEC C***CATEGORY K1A1A1, K1A2A, L8A3 C***TYPE DOUBLE PRECISION (FC-S, DFC-D) C***KEYWORDS B-SPLINE, CONSTRAINED LEAST SQUARES, CURVE FITTING, C 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 Equality and inequality constraints can be imposed on the C fitted curve. 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 documentation and discussion of constrained C curve fitting using 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 DFC( ) 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 NCONST,XCONST(*), C YCONST(*),NDERIV(*) C The number of conditions that constrain the C B-spline is NCONST. A constraint is specified C by an (X,Y) pair in the arrays XCONST(*) and C YCONST(*), and by the type of constraint and C derivative value encoded in the array C NDERIV(*). No sorting of XCONST(*) is C required. The value of NDERIV(*) is C determined as follows. Suppose the I-th C constraint applies to the J-th derivative C of the B-spline. (Any non-negative value of C J < NORD is permitted. In particular the C value J=0 refers to the B-spline itself.) C For this I-th constraint, set C XCONST(I)=X, C YCONST(I)=Y, and C NDERIV(I)=ITYPE+4*J, where C C ITYPE = 0, if (J-th deriv. at X) .LE. Y. C = 1, if (J-th deriv. at X) .GE. Y. C = 2, if (J-th deriv. at X) .EQ. Y. C = 3, if (J-th deriv. at X) .EQ. C (J-th deriv. at Y). C (A value of NDERIV(I)=-1 will cause this C constraint to be ignored. This subprogram C feature is often useful when temporarily C suppressing a constraint while still C retaining the source code of the calling C program.) C C MODE C An input flag that directs the least squares C solution method used by DFC( ). C C The variance function, referred to below, C defines the square of the probable error of C the fitted curve at any point, XVAL. C This feature of DFC( ) allows one to use the C square root of this variance function to C determine a probable error band around the C fitted curve. C C =1 a new problem. No variance function. C C =2 a new problem. Want variance function. C C =3 an old problem. No variance function. C C =4 an old problem. Want variance function. C C Any value of MODE other than 1-4 is an error. C C The user with a new problem can skip directly C to the description of the input parameters C IW(1), IW(2). C C If the user correctly specifies the new or old C problem status, the subprogram DFC( ) will C perform more efficiently. C By an old problem it is meant that subprogram C DFC( ) was last called with this same set of C knots, data points and weights. C C Another often useful deployment of this old C problem designation can occur when one has C previously obtained a Q-R orthogonal C decomposition of the matrix resulting from C B-spline fitting of data (without constraints) C at the breakpoints BKPT(I), I=1,...,NBKPT. C For example, this matrix could be the result C of sequential accumulation of the least C squares equations for a very large data set. C The user writes this code in a manner C convenient for the application. For the C discussion here let C C N=NBKPT-NORD, and K=N+3 C C Let us assume that an equivalent least squares C system C C RC=D C C has been obtained. Here R is an N+1 by N C matrix and D is a vector with N+1 components. C The last row of R is zero. The matrix R is C upper triangular and banded. At most NORD of C the diagonals are nonzero. C The contents of R and D can be copied to the C working array W(*) as follows. C C The I-th diagonal of R, which has N-I+1 C elements, is copied to W(*) starting at C C W((I-1)*K+1), C C for I=1,...,NORD. C The vector D is copied to W(*) starting at C C W(NORD*K+1) C C The input value used for NDATA is arbitrary C when an old problem is designated. Because C of the feature of DFC( ) that checks the C working storage array lengths, a value not C exceeding NBKPT should be used. For example, C use NDATA=0. C C (The constraints or variance function request C can change in each call to DFC( ).) A new C problem is anything other than an old problem. C C IW(1),IW(2) C The amounts of working storage actually C allocated for the working arrays W(*) and C IW(*). These quantities are compared with the C actual amounts of storage needed in DFC( ). C Insufficient storage allocated for either C W(*) or IW(*) is an error. This feature was C included in DFC( ) because misreading the C storage formulas for W(*) and IW(*) might very C well lead to subtle and hard-to-find C programming bugs. C C The length of W(*) must be at least C C NB=(NBKPT-NORD+3)*(NORD+1)+ C 2*MAX(NDATA,NBKPT)+NBKPT+NORD**2 C C Whenever possible the code uses banded matrix C processors DBNDAC( ) and DBNDSL( ). These C are utilized if there are no constraints, C no variance function is required, and there C is sufficient data to uniquely determine the C B-spline coefficients. If the band processors C cannot be used to determine the solution, C then the constrained least squares code DLSEI C is used. In this case the subprogram requires C an additional block of storage in W(*). For C the discussion here define the integers NEQCON C and NINCON respectively as the number of C equality (ITYPE=2,3) and inequality C (ITYPE=0,1) constraints imposed on the fitted C curve. Define C C L=NBKPT-NORD+1 C C and note that C C NCONST=NEQCON+NINCON. C C When the subprogram DFC( ) uses DLSEI( ) the C length of the working array W(*) must be at C least C C LW=NB+(L+NCONST)*L+ C 2*(NEQCON+L)+(NINCON+L)+(NINCON+2)*(L+6) C C The length of the array IW(*) must be at least C C IW1=NINCON+2*L C C in any case. C C Output.. All TYPE REAL variables are DOUBLE PRECISION C MODE C An output flag that indicates the status C of the constrained curve fit. C C =-1 a usage error of DFC( ) occurred. The C offending condition is noted with the C SLATEC library error processor, XERMSG. C In case the working arrays W(*) or IW(*) C are not long enough, the minimal C acceptable length is printed. C C = 0 successful constrained curve fit. C C = 1 the requested equality constraints C are contradictory. C C = 2 the requested inequality constraints C are contradictory. C C = 3 both equality and inequality constraints C are contradictory. C C COEFF(*) C If the output value of MODE=0 or 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 MODE=1, the equality constraints are C contradictory. To make the fitting process C more robust, the equality constraints are C satisfied in a least squares sense. In this C case the array COEFF(*) contains B-spline C coefficients for this extended concept of a C solution. If MODE=-1,2 or 3 on output, the C array COEFF(*) is undefined. C C Working Arrays.. All Type REAL variables are DOUBLE PRECISION C W(*),IW(*) C These arrays are respectively typed DOUBLE C PRECISION and INTEGER. C Their required lengths are specified as input C parameters in IW(1), IW(2) noted above. The C contents of W(*) must not be modified by the C user if the variance function is desired. C C Evaluating the C Variance Function.. C To evaluate the variance function (assuming C that the uncertainties of the Y values were C provided to DFC( ) and an input value of C MODE=2 or 4 was used), use the function C subprogram DCV( ) C C VAR=DCV(XVAL,NDATA,NCONST,NORD,NBKPT, C BKPT,W) C C Here XVAL is the point where the variance is C desired. The other arguments have the same C meaning as in the usage of DFC( ). C C For those users employing the old problem C designation, let MDATA be the number of data C points in the problem. (This may be different C from NDATA if the old problem designation C feature was used.) The value, VAR, should be C multiplied by the quantity C C DBLE(MAX(NDATA-N,1))/DBLE(MAX(MDATA-N,1)) C C The output of this subprogram is not defined C if an input value of MODE=1 or 3 was used in C FC( ) or if an output value of MODE=-1, 2, or C 3 was obtained. The variance function, except C for the scaling factor noted above, is given C by C C VAR=(transpose of B(XVAL))*C*B(XVAL) C C The vector B(XVAL) is the B-spline basis C function values at X=XVAL. C The covariance matrix, C, of the solution C coefficients accounts only for the least C squares equations and the explicitly stated C equality constraints. This fact must be C considered when interpreting the variance C function from a data fitting problem that has C inequality constraints on the fitted curve. 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 MODE=0 or 1 C was obtained from DFC( ), 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 C the call to DFC( ). 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 DFCMN C***REVISION HISTORY (YYMMDD) C 780801 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 891006 Cosmetic changes to prologue. (WRB) C 891006 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900510 Convert references to XERRWV to references 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 DFC DOUBLE PRECISION BKPT(*), COEFF(*), SDDATA(*), W(*), XCONST(*), * XDATA(*), YCONST(*), YDATA(*) INTEGER IW(*), MODE, NBKPT, NCONST, NDATA, NDERIV(*), NORD C EXTERNAL DFCMN C INTEGER I1, I2, I3, I4, I5, I6, I7, MDG, MDW C C***FIRST EXECUTABLE STATEMENT DFC MDG = NBKPT - NORD + 3 MDW = NBKPT - NORD + 1 + NCONST C USAGE IN DFCMN( ) OF W(*).. C I1,...,I2-1 G(*,*) C C I2,...,I3-1 XTEMP(*) C C I3,...,I4-1 PTEMP(*) C C I4,...,I5-1 BKPT(*) (LOCAL TO DFCMN( )) C C I5,...,I6-1 BF(*,*) C C I6,...,I7-1 W(*,*) C C I7,... WORK(*) FOR DLSEI( ) C I1 = 1 I2 = I1 + MDG*(NORD+1) I3 = I2 + MAX(NDATA,NBKPT) I4 = I3 + MAX(NDATA,NBKPT) I5 = I4 + NBKPT I6 = I5 + NORD*NORD I7 = I6 + MDW*(NBKPT-NORD+1) CALL DFCMN(NDATA, XDATA, YDATA, SDDATA, NORD, NBKPT, BKPT, NCONST, 1 XCONST, YCONST, NDERIV, MODE, COEFF, W(I5), W(I2), W(I3), 2 W(I4), W(I1), MDG, W(I6), MDW, W(I7), IW) RETURN END