SUBROUTINE FC(NDATA,XDATA,YDATA,SDDATA,NORD,NBKPT,BKPT,NCONST, 1 XCONST,YCONST,NDERIV,MODE,COEFF,W,IW) C***BEGIN PROLOGUE FC C***DATE WRITTEN 780801 (YYMMDD) C***REVISION DATE 820801 (YYMMDD) C***CATEGORY NO. K1A1A1,K1A2A,L8A3 C***KEYWORDS B-SPLINES,CONSTRAINED LEAST SQUARES,CURVE FITTING, C LEAST SQUARES C***AUTHOR HANSON, R. J., (SNLA) C***PURPOSE Fits a piece-wise polynomial curve to discrete data. C The piece-wise polynomials are represented as B-splines. C The fitting is done in a least squares sense. Equality C and inequality constraints can be imposed on the fitted C curve. C***DESCRIPTION C C This subprogram fits a piece-wise polynomial curve C to discrete data. The piece-wise 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. For a description of the B-splines and C usage instructions and software for evaluating 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.. C NDATA,XDATA(*), C YDATA(*), C SDDATA(*) C The NDATA discrete (X,Y) pairs and C the Y value standard deviation or C uncertainty, SD, are in the respective C arrays XDATA(*), YDATA(*), and SDDATA(*). C No sorting of XDATA(*) is required. Any C non-negative value of NDATA is allowed. A C negative value of NDATA is an error. C A zero value for any entry of SDDATA(*) C 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 C NORD are in the array BKPT(*). Normally C the problem data interval will be included C between the limits BKPT(NORD) and C BKPT(NBKPT-NORD+1). The additional end knots C BKPT(I),I=1,...,NORD-1 and I=NBKPT-NORD+2,..., C NBKPT, are required to compute the functions C used to fit the data. C No sorting of BKPT(*) is required. C Internal to FC( ) the extreme end knots C may be slightly reduced and increased C respectively to accommodate any data C values that are exterior to the given knot C values. The array BKPT(*) is not altered. C The value of NORD must satisfy 1 < NORD < 20 . C The value of NBKPT must satisfy NBKPT .GE. C 2*NORD. Other values are considered errors. C C (The order of the spline is one more C than the degree of the piece-wise polynomial C defined on each interval. This is consistent C with the B-spline package convention. For C example, NORD=4 when we are using piece-wise C cubics.) C C NCONST,XCONST(*), C YCONST(*),NDERIV(*) C The number of conditions that constrain C the B-spline is NCONST. A constraint is C specified by an (X,Y) pair in the arrays C XCONST(*) and YCONST(*), and by the type C of constraint and derivative value encoded C in the array NDERIV(*). No sorting C of XCONST(*) is required. The value of C NDERIV(*) is determined as follows. C Suppose the I-th constraint applies to the C J-th derivative of the B-spline. (Any C non-negative value of J < NORD is permitted. C In particular the value J=0 refers to the C 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 C squares solution method used by FC( ). C C The variance function, referred to below, C defines the square of the probable error C of the fitted curve at any point, XVAL. C This feature of FC( ) allows one to C use the square root of this variance C function to determine a probable error C band around the 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 C directly to the description of the C input parameters IW(1), IW(2). C C C If the user correctly specifies the new or old C problem status, the subprogram FC( ) will C perform more efficiently. C By an old problem it is meant that subprogram C FC( ) was last called with this same set C of knots, data points and weights. C C Another often useful deployment of this C old problem designation can occur when one C has previously obtained a Q-R orthogonal C decomposition of the matrix resulting C from B-spline fitting of data (without C constraints) at the breakpoints BKPT(I), C I=1,...,NBKPT. For example, this matrix C could be the result of sequential C accumulation of the least squares C 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 C squares system C C RC=D C C has been obtained. Here R is an N+1 C by N matrix and D is a vector with C N+1 components. The last row of R C is zero. The matrix R is upper triangular C and banded. At most NORD of the diagonals are C nonzero. C The contents of R and D can be copied to C the working array W(*) as follows. C C The I-th diagonal of R, which has C N-I+1 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 FC( ) that checks the C working storage array lengths, a value C not exceeding NBKPT should be used. C For example, use NDATA=0. C C (The constraints or variance function C request can change in each call to FC( ).) C A new problem is anything other than an old C problem. C 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 FC( ). C Insufficient storage allocated for C either W(*) or IW(*) is an error. C This feature was included in FC( ) C because misreading the storage formulas C for W(*) and IW(*) might very well C lead to subtle and hard-to-find C programming bugs. C C Length of W(*) must be at least C C NB=(NBKPT-NORD+3)*(NORD+1)+ C 2*MAX0(NDATA,NBKPT)+NBKPT+NORD**2 C C Whenever possible the code uses banded matrix C processors BNDACC( ) and BNDSOL( ). These C are utilized if there are no constraints, C no variance function is required, and there C is sufficient data to uniquely determine C the B-spline coefficents. C If the band processors cannot be used to C determine the solution, then the constrained C least squares code LSEI( ) is used. C In this case the subprogram requires an C additional block of storage in W(*). For the C discussion here define the integers C NEQCON and NINCON respectively as the C number of equality (ITYPE=2,3) and C inequality (ITYPE=0,1) constraints C imposed on the fitted curve. Define C C L=NBKPT-NORD+1 C C and note that C C NCONST=NEQCON+NINCON. C C When the subprogram FC( ) uses LSEI( ) 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 Output.. C MODE C An output flag that indicates the status C of the constrained curve fit. C C =-1 a usage error of FC( ) occurred. The C offending condition is noted with the C SLATEC library error processor, C XERROR( ). In case the working C arrays W(*) or IW(*) are not long C enough, the minimal acceptable C length is printed using the SLATEC C error processing subprogram, XERRWV( ). 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 C this case the array COEFF(*) contains C B-spline coefficients for this extended C concept of a solution. C If MODE=-1,2 or 3 on output, the array C COEFF(*) is undefined. C C Working Arrays.. C W(*),IW(*) C These arrays are respectively typed C real and integer. Their required C lengths are specified as input parameters C in IW(1), IW(2) noted above. C The contents of W(*) must not be modified C by the user if the variance function C is desired. C C Evaluating the C Variance Function.. C To evaluate the variance function C (assuming that the uncertainties of the C Y values were provided to FC( ) and an C input value of MODE=2 or 4 was used), use the C function subprogram CV( ) C C VAR=CV(XVAL,NDATA,NCONST,NORD,NBKPT, C BKPT,W) C C Here XVAL is the point where the variance C is desired. The other arguments have the C same meaning as in the usage of FC( ). C C For those users employing the old problem C designation, let MDATA be the number of C data points in the problem. (This C may be different from NDATA if the C old problem designation feature was C used.) The value, VAR, should be multiplied C by the quantity C C FLOAT(MAX0(NDATA-N,1))/FLOAT(MAX0(MDATA-N,1)) C C The output of this subprogram is not C defined if an input value of MODE=1 or 3 C was used in FC( ) or if an output value of C MODE=-1, 2, or 3 was obtained. C The variance function, except for the C scaling factor noted above, is given 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 C solution coefficients accounts only C for the least squares equations and C the explicitly stated equality constraints. C This fact must be considered when C interpreting the variance function from C 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 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 MODE=0 or 1 C was obtained from FC( ), 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 C the call to FC( ). C C BVALU( ) expects the breakpoint array BKPT(*) C to be sorted. C***REFERENCES HANSON R.J., *CONSTRAINED LEAST SQUARES CURVE FITTING C TO DISCRETE DATA USING B-SPLINES, A USERS C GUIDE*, SAND78-1291, DECEMBER,1978. C***ROUTINES CALLED FCMN C***END PROLOGUE FC