SLATEC Routines --- DFC ---
*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.class" tppabs="http://www.kachinatech.com/~hjjou/slatec-doc/code/of.class" 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