This subprogram fits a piecewise polynomial curve to discrete data. The piecewise polynomials are represented as Bsplines. The fitting is done in a weighted least squares sense. Equality and inequality constraints can be imposed on the fitted curve.
To evaluate the variance function (assuming that the uncertainties of the Y values were provided to DFC and an input value of MODE=2 or 4 was used), use the function subprogram DCV
var = dcv(xval,ndata,nconst,nord,nbkpt, bkpt,w)
Here XVAL is the point where the variance is desired. The other arguments have the same meaning as in the usage of DFC.
For those users employing the old problem designation, let MDATA be the number of data points in the problem. (This may be different from NDATA if the old problem designation feature was used.) The value, VAR, should be multiplied by the quantity
DBLE(MAX(NDATAN,1))/DBLE(MAX(MDATAN,1))
The output of this subprogram is not defined if an input value of MODE=1 or 3 was used in FC( ) or if an output value of MODE=1, 2, or 3 was obtained. The variance function, except for the scaling factor noted above, is given by
VAR=(transpose of B(XVAL))*C*B(XVAL)
The vector B(XVAL) is the Bspline basis function values at X=XVAL. The covariance matrix, C, of the solution coefficients accounts only for the least squares equations and the explicitly stated equality constraints. This fact must be considered when interpreting the variance function from a data fitting problem that has inequality constraints on the fitted curve.
Type  Intent  Optional  Attributes  Name  

integer,  intent(in)  ::  ndata 
number of points (size of 

real(kind=wp),  intent(in)  ::  xdata(*) 
X data array. No sorting of 

real(kind=wp),  intent(in)  ::  ydata(*) 
Y data array. 

real(kind=wp),  intent(in)  ::  sddata(*) 
Y value standard deviation or uncertainty.
A zero value for any entry of


integer,  intent(in)  ::  nord 
Bspline order.
(The order of the spline is one more than the
degree of the piecewise polynomial defined on
each interval. This is consistent with the
Bspline package convention. For example,


integer,  intent(in)  ::  Nbkpt 
The value of 

real(kind=wp),  intent(in),  dimension(*)  ::  Bkpt 


integer,  intent(in)  ::  nconst 
The number of conditions that constrain the Bspline is NCONST. A constraint is specified by an (X,Y) pair in the arrays XCONST() and YCONST(), and by the type of constraint and derivative value encoded in the array NDERIV(*). 

real(kind=wp),  intent(in)  ::  xconst(*) 
X value of constraint. No sorting of XCONST(*) is required. 

real(kind=wp),  intent(in)  ::  yconst(*) 
Y value of constraint 

integer,  intent(in)  ::  nderiv(*) 
The value of NDERIV(*) is determined as follows. Suppose the Ith constraint applies to the Jth derivative of the Bspline. (Any nonnegative value of J < NORD is permitted. In particular the value J=0 refers to the Bspline itself.) For this Ith constraint, set
(A value of NDERIV(I)=1 will cause this constraint to be ignored. This subprogram feature is often useful when temporarily suppressing a constraint while still retaining the source code of the calling program.) 

integer,  intent(inout)  ::  mode 
Input An input flag that directs the least squares solution method used by DFC. The variance function, referred to below, defines the square of the probable error of the fitted curve at any point, XVAL. This feature of DFC allows one to use the square root of this variance function to determine a probable error band around the fitted curve.
Any value of MODE other than 14 is an error. The user with a new problem can skip directly to the description of the input parameters IW(1), IW(2). If the user correctly specifies the new or old problem status, the subprogram DFC will perform more efficiently. By an old problem it is meant that subprogram DFC was last called with this same set of knots, data points and weights. Another often useful deployment of this old problem designation can occur when one has previously obtained a QR orthogonal decomposition of the matrix resulting from Bspline fitting of data (without constraints) at the breakpoints BKPT(I), I=1,...,NBKPT. For example, this matrix could be the result of sequential accumulation of the least squares equations for a very large data set. The user writes this code in a manner convenient for the application. For the discussion here let
Let us assume that an equivalent least squares system
has been obtained. Here R is an N+1 by N matrix and D is a vector with N+1 components. The last row of R is zero. The matrix R is upper triangular and banded. At most NORD of the diagonals are nonzero. The contents of R and D can be copied to the working array W(*) as follows. The Ith diagonal of R, which has NI+1 elements, is copied to W(*) starting at
for I=1,...,NORD. The vector D is copied to W(*) starting at
The input value used for NDATA is arbitrary when an old problem is designated. Because of the feature of DFC that checks the working storage array lengths, a value not exceeding NBKPT should be used. For example, use NDATA=0. (The constraints or variance function request can change in each call to DFC.) A new problem is anything other than an old problem. Output An output flag that indicates the status of the constrained curve fit.


real(kind=wp),  intent(out)  ::  coeff(*) 
If the output value of MODE=0 or 1, this array contains the unknowns obtained from the least squares fitting process. These N=NBKPTNORD parameters are the Bspline coefficients. For MODE=1, the equality constraints are contradictory. To make the fitting process more robust, the equality constraints are satisfied in a least squares sense. In this case the array COEFF() contains Bspline coefficients for this extended concept of a solution. If MODE=1,2 or 3 on output, the array COEFF() is undefined. 

real(kind=wp)  ::  w(*) 
real work array of length The length of W(*) must be at least
Whenever possible the code uses banded matrix processors DBNDAC( ) and DBNDSL( ). These are utilized if there are no constraints, no variance function is required, and there is sufficient data to uniquely determine the Bspline coefficients. If the band processors cannot be used to determine the solution, then the constrained least squares code DLSEI is used. In this case the subprogram requires an additional block of storage in W(*). For the discussion here define the integers NEQCON and NINCON respectively as the number of equality (ITYPE=2,3) and inequality (ITYPE=0,1) constraints imposed on the fitted curve. Define
and note that
When the subprogram DFC uses DLSEI the length of the working array W(*) must be at least


integer  ::  iw(*) 
integer work array of length
The length of the array IW(*) must be at least
in any case. 
subroutine dfc (ndata, xdata, ydata, sddata, nord, nbkpt, bkpt, & nconst, xconst, yconst, nderiv, mode, coeff, w, iw) integer, intent(in) :: ndata !! number of points (size of `xdata` and `ydata`). !! Any nonnegative value of `NDATA` is allowed. !! A negative value of `NDATA` is an error. real(wp), intent(in) :: xdata(*) !! X data array. No sorting of `XDATA(*)` is required. real(wp), intent(in) :: ydata(*) !! Y data array. real(wp), intent(in) :: sddata(*) !! Y value standard deviation or uncertainty. !! A zero value for any entry of !! `SDDATA(*)` will weight that data point as 1. !! Otherwise the weight of that data point is !! the reciprocal of this entry. integer, intent(in) :: nord !! Bspline order. !! (The order of the spline is one more than the !! degree of the piecewise polynomial defined on !! each interval. This is consistent with the !! Bspline package convention. For example, !! `NORD=4` when we are using piecewise cubics.) !! `NORD` must be in the range `1 <= NORD <= 20`. integer,intent(in) :: Nbkpt !! The value of `NBKPT` must satisfy the condition `NBKPT >= 2*NORD`. real(wp),dimension(*),intent(in) :: Bkpt !! `NBKPT` knots of the Bspline. !! Normally the !! problem data interval will be included between !! the limits `BKPT(NORD)` and `BKPT(NBKPTNORD+1)`. !! The additional end knots `BKPT(I),I=1,...,NORD1` !! and `I=NBKPTNORD+2,...,NBKPT`, are !! required to compute the functions used to fit !! the data. No sorting of `BKPT(*)` is required. !! Internal to [[DEFC]] the extreme end knots may !! be reduced and increased respectively to !! accommodate any data values that are exterior !! to the given knot values. The contents of !! `BKPT(*)` is not changed. integer, intent(in) :: nconst !! The number of conditions that constrain the !! Bspline is NCONST. A constraint is specified !! by an (X,Y) pair in the arrays XCONST(*) and !! YCONST(*), and by the type of constraint and !! derivative value encoded in the array !! NDERIV(*). real(wp), intent(in) :: xconst(*) !! X value of constraint. !! No sorting of XCONST(*) is required. real(wp), intent(in) :: yconst(*) !! Y value of constraint integer, intent(in) :: nderiv(*) !! The value of NDERIV(*) is !! determined as follows. Suppose the Ith !! constraint applies to the Jth derivative !! of the Bspline. (Any nonnegative value of !! J < NORD is permitted. In particular the !! value J=0 refers to the Bspline itself.) !! For this Ith constraint, set !!``` !! XCONST(I)=X, !! YCONST(I)=Y, and !! NDERIV(I)=ITYPE+4*J, where !! !! ITYPE = 0, if (Jth deriv. at X) <= Y. !! = 1, if (Jth deriv. at X) >= Y. !! = 2, if (Jth deriv. at X) == Y. !! = 3, if (Jth deriv. at X) == !! (Jth deriv. at Y). !!``` !! (A value of NDERIV(I)=1 will cause this !! constraint to be ignored. This subprogram !! feature is often useful when temporarily !! suppressing a constraint while still !! retaining the source code of the calling !! program.) integer, intent(inout) :: mode !! *Input* !! !! An input flag that directs the least squares !! solution method used by [[DFC]]. !! !! The variance function, referred to below, !! defines the square of the probable error of !! the fitted curve at any point, XVAL. !! This feature of [[DFC]] allows one to use the !! square root of this variance function to !! determine a probable error band around the !! fitted curve. !! !! * `=1` a new problem. No variance function. !! * `=2` a new problem. Want variance function. !! * `=3` an old problem. No variance function. !! * `=4` an old problem. Want variance function. !! !! Any value of MODE other than 14 is an error. !! !! The user with a new problem can skip directly !! to the description of the input parameters !! IW(1), IW(2). !! !! If the user correctly specifies the new or old !! problem status, the subprogram [[DFC]] will !! perform more efficiently. !! By an old problem it is meant that subprogram !! [[DFC]] was last called with this same set of !! knots, data points and weights. !! !! Another often useful deployment of this old !! problem designation can occur when one has !! previously obtained a QR orthogonal !! decomposition of the matrix resulting from !! Bspline fitting of data (without constraints) !! at the breakpoints BKPT(I), I=1,...,NBKPT. !! For example, this matrix could be the result !! of sequential accumulation of the least !! squares equations for a very large data set. !! The user writes this code in a manner !! convenient for the application. For the !! discussion here let !! !! `N=NBKPTNORD, and K=N+3` !! !! Let us assume that an equivalent least squares !! system !! !! `RC=D` !! !! has been obtained. Here R is an N+1 by N !! matrix and D is a vector with N+1 components. !! The last row of R is zero. The matrix R is !! upper triangular and banded. At most NORD of !! the diagonals are nonzero. !! The contents of R and D can be copied to the !! working array W(*) as follows. !! !! The Ith diagonal of R, which has NI+1 !! elements, is copied to W(*) starting at !! !! `W((I1)*K+1),` !! !! for I=1,...,NORD. !! The vector D is copied to W(*) starting at !! !! `W(NORD*K+1)` !! !! The input value used for NDATA is arbitrary !! when an old problem is designated. Because !! of the feature of [[DFC]] that checks the !! working storage array lengths, a value not !! exceeding NBKPT should be used. For example, !! use NDATA=0. !! !! (The constraints or variance function request !! can change in each call to [[DFC]].) A new !! problem is anything other than an old problem. !! !! *Output* !! !! An output flag that indicates the status !! of the constrained curve fit. !! !! * `=1` a usage error of [[DFC]] occurred. The !! offending condition is noted with the !! SLATEC library error processor, XERMSG. !! In case the working arrays W(*) or IW(*) !! are not long enough, the minimal !! acceptable length is printed. !! * `= 0` successful constrained curve fit. !! * `= 1` the requested equality constraints !! are contradictory. !! * `= 2` the requested inequality constraints !! are contradictory. !! * `= 3` both equality and inequality constraints !! are contradictory. real(wp), intent(out) :: coeff(*) !! If the output value of MODE=0 or 1, this array !! contains the unknowns obtained from the least !! squares fitting process. These N=NBKPTNORD !! parameters are the Bspline coefficients. !! For MODE=1, the equality constraints are !! contradictory. To make the fitting process !! more robust, the equality constraints are !! satisfied in a least squares sense. In this !! case the array COEFF(*) contains Bspline !! coefficients for this extended concept of a !! solution. If MODE=1,2 or 3 on output, the !! array COEFF(*) is undefined. real(wp) :: w(*) !! real work array of length `IW(1)`. The !! contents of `W(*)` must not be modified by the !! user if the variance function is desired. !! !! The length of W(*) must be at least !!``` !! NB=(NBKPTNORD+3)*(NORD+1)+ !! 2*MAX(NDATA,NBKPT)+NBKPT+NORD**2 !!``` !! Whenever possible the code uses banded matrix !! processors DBNDAC( ) and DBNDSL( ). These !! are utilized if there are no constraints, !! no variance function is required, and there !! is sufficient data to uniquely determine the !! Bspline coefficients. If the band processors !! cannot be used to determine the solution, !! then the constrained least squares code DLSEI !! is used. In this case the subprogram requires !! an additional block of storage in W(*). For !! the discussion here define the integers NEQCON !! and NINCON respectively as the number of !! equality (ITYPE=2,3) and inequality !! (ITYPE=0,1) constraints imposed on the fitted !! curve. Define !! !! `L = NBKPTNORD+1` !! !! and note that !! !! `NCONST = NEQCON+NINCON` !! !! When the subprogram [[DFC]] uses [[DLSEI]] the !! length of the working array W(*) must be at !! least !! !! `LW = NB+(L+NCONST)*L+2*(NEQCON+L)+(NINCON+L)+(NINCON+2)*(L+6)` integer :: iw(*) !! integer work array of length `IW(2)` !! !! `IW(1),IW(2)` are the amounts of working storage actually !! allocated for the working arrays W(*) and !! IW(*). These quantities are compared with the !! actual amounts of storage needed in [[DFC]]. !! Insufficient storage allocated for either !! W(*) or IW(*) is an error. This feature was !! included in [[DFC]] because misreading the !! storage formulas for W(*) and IW(*) might very !! well lead to subtle and hardtofind !! programming bugs. !! !! The length of the array IW(*) must be at least !! !! `IW1 = NINCON+2*L` !! !! in any case. integer :: i1, i2, i3, i4, i5, i6, i7, mdg, mdw mdg = nbkpt  nord + 3 mdw = nbkpt  nord + 1 + nconst ! USAGE IN DFCMN( ) OF W(*).. ! I1,...,I21 G(*,*) ! I2,...,I31 XTEMP(*) ! I3,...,I41 PTEMP(*) ! I4,...,I51 BKPT(*) (LOCAL TO [[DFCMN]]) ! I5,...,I61 BF(*,*) ! I6,...,I71 W(*,*) ! I7,... WORK(*) FOR [[DLSEI]] 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*(nbkptnord+1) call dfcmn(ndata, xdata, ydata, sddata, nord, nbkpt, bkpt, nconst, & xconst, yconst, nderiv, mode, coeff, w(i5), w(i2), w(i3), & w(i4), w(i1), mdg, w(i6), mdw, w(i7), iw) end subroutine dfc