dfc Subroutine

public subroutine dfc(ndata, xdata, ydata, sddata, nord, Nbkpt, Bkpt, nconst, xconst, yconst, nderiv, mode, coeff, w, iw)

This subprogram fits a piecewise polynomial curve to discrete data. The piecewise polynomials are represented as B-splines. The fitting is done in a weighted least squares sense. Equality and inequality constraints can be imposed on the fitted curve.

Evaluating the Variance Function

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(NDATA-N,1))/DBLE(MAX(MDATA-N,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 B-spline 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.

Evaluating the Fitted Curve

  • Refer to the defc header

Revision history

  • 780801 DATE WRITTEN. Hanson, R. J., (SNLA)
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 891006 Cosmetic changes to prologue. (WRB)
  • 891006 REVISION DATE from Version 3.2
  • 891214 Prologue converted to Version 4.0 format. (BAB)
  • 900510 Convert references to XERRWV to references to XERMSG. (RWC)
  • 900607 Editorial changes to Prologue to make Prologues for EFC, DEFC, FC, and DFC look as much the same as possible. (RWC)
  • 920501 Reformatted the REFERENCES section. (WRB)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: ndata

number of points (size of xdata and ydata). Any non-negative value of NDATA is allowed. A negative value of NDATA is an error.
real(kind=wp), intent(in) :: xdata(*)

X data array. No sorting of XDATA(*) is required.
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 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

B-spline 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 B-spline 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(kind=wp), intent(in), dimension(*) :: Bkpt

NBKPT knots of the B-spline. Normally the problem data interval will be included between the limits BKPT(NORD) and BKPT(NBKPT-NORD+1). The additional end knots BKPT(I),I=1,...,NORD-1 and I=NBKPT-NORD+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 B-spline 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 I-th constraint applies to the J-th derivative of the B-spline. (Any non-negative value of J < NORD is permitted. In particular the value J=0 refers to the B-spline itself.) For this I-th constraint, set

  XCONST(I)=X,
  YCONST(I)=Y, and
  NDERIV(I)=ITYPE+4*J, where

  ITYPE = 0,      if (J-th deriv. at X) <= Y.
        = 1,      if (J-th deriv. at X) >= Y.
        = 2,      if (J-th deriv. at X) == Y.
        = 3,      if (J-th deriv. at X) ==
                     (J-th 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 1-4 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 Q-R orthogonal decomposition of the matrix resulting from B-spline 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=NBKPT-NORD, 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 I-th diagonal of R, which has N-I+1 elements, is copied to W(*) starting at

W((I-1)*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(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=NBKPT-NORD parameters are the B-spline 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 B-spline 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 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=(NBKPT-NORD+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 B-spline 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 = NBKPT-NORD+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 hard-to-find programming bugs.

The length of the array IW(*) must be at least

IW1 = NINCON+2*L

in any case.

Calls

proc~~dfc~~CallsGraph proc~dfc bspline_defc_module::dfc proc~dfcmn bspline_defc_module::dfcmn proc~dfc->proc~dfcmn proc~daxpy bspline_blas_module::daxpy proc~dfcmn->proc~daxpy proc~dbndac bspline_defc_module::dbndac proc~dfcmn->proc~dbndac proc~dbndsl bspline_defc_module::dbndsl proc~dfcmn->proc~dbndsl proc~dcopy bspline_blas_module::dcopy proc~dfcmn->proc~dcopy proc~dfspvd bspline_defc_module::dfspvd proc~dfcmn->proc~dfspvd proc~dfspvn bspline_defc_module::dfspvn proc~dfcmn->proc~dfspvn proc~dlsei bspline_defc_module::dlsei proc~dfcmn->proc~dlsei proc~dscal bspline_blas_module::dscal proc~dfcmn->proc~dscal proc~dsort bspline_defc_module::dsort proc~dfcmn->proc~dsort proc~dh12 bspline_defc_module::dh12 proc~dbndac->proc~dh12 proc~dfspvd->proc~dfspvn proc~dlsei->proc~daxpy proc~dlsei->proc~dcopy proc~dlsei->proc~dscal proc~dasum bspline_blas_module::dasum proc~dlsei->proc~dasum proc~ddot bspline_blas_module::ddot proc~dlsei->proc~ddot proc~dlsei->proc~dh12 proc~dlsi bspline_defc_module::dlsi proc~dlsei->proc~dlsi proc~dnrm2 bspline_blas_module::dnrm2 proc~dlsei->proc~dnrm2 proc~dswap bspline_blas_module::dswap proc~dlsei->proc~dswap proc~sort_ascending bspline_defc_module::sort_ascending proc~dsort->proc~sort_ascending proc~dh12->proc~daxpy proc~dh12->proc~ddot proc~dh12->proc~dswap proc~dlsi->proc~daxpy proc~dlsi->proc~dcopy proc~dlsi->proc~dscal proc~dlsi->proc~dasum proc~dlsi->proc~ddot proc~dlsi->proc~dh12 proc~dlsi->proc~dswap proc~dhfti bspline_defc_module::dhfti proc~dlsi->proc~dhfti proc~dlpdp bspline_defc_module::dlpdp proc~dlsi->proc~dlpdp proc~dhfti->proc~dh12 proc~dlpdp->proc~dcopy proc~dlpdp->proc~dscal proc~dlpdp->proc~ddot proc~dlpdp->proc~dnrm2 proc~dwnnls bspline_defc_module::dwnnls proc~dlpdp->proc~dwnnls proc~dwnlsm bspline_defc_module::dwnlsm proc~dwnnls->proc~dwnlsm proc~dwnlsm->proc~daxpy proc~dwnlsm->proc~dcopy proc~dwnlsm->proc~dscal proc~dwnlsm->proc~dasum proc~dwnlsm->proc~dh12 proc~dwnlsm->proc~dnrm2 proc~dwnlsm->proc~dswap proc~drotm bspline_blas_module::drotm proc~dwnlsm->proc~drotm proc~drotmg bspline_blas_module::drotmg proc~dwnlsm->proc~drotmg proc~dwnlit bspline_defc_module::dwnlit proc~dwnlsm->proc~dwnlit proc~idamax bspline_blas_module::idamax proc~dwnlsm->proc~idamax proc~dwnlit->proc~dcopy proc~dwnlit->proc~dscal proc~dwnlit->proc~dh12 proc~dwnlit->proc~dswap proc~dwnlit->proc~drotm proc~dwnlit->proc~drotmg proc~dwnlit->proc~idamax proc~dwnlt1 bspline_defc_module::dwnlt1 proc~dwnlit->proc~dwnlt1 proc~dwnlt2 bspline_defc_module::dwnlt2 proc~dwnlit->proc~dwnlt2 proc~dwnlt3 bspline_defc_module::dwnlt3 proc~dwnlit->proc~dwnlt3 proc~dwnlt1->proc~idamax proc~dwnlt3->proc~dswap

Source Code

   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 non-negative 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 !! B-spline 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
                               !! B-spline 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 B-spline.
                                            !! Normally the
                                            !! problem data interval will be included between
                                            !! the limits `BKPT(NORD)` and `BKPT(NBKPT-NORD+1)`.
                                            !! The additional end knots `BKPT(I),I=1,...,NORD-1`
                                            !! and `I=NBKPT-NORD+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
                                 !! B-spline 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 I-th
                                    !! constraint applies to the J-th derivative
                                    !! of the B-spline.  (Any non-negative value of
                                    !! J < NORD is permitted.  In particular the
                                    !! value J=0 refers to the B-spline itself.)
                                    !! For this I-th constraint, set
                                    !!```
                                    !!  XCONST(I)=X,
                                    !!  YCONST(I)=Y, and
                                    !!  NDERIV(I)=ITYPE+4*J, where
                                    !!
                                    !!  ITYPE = 0,      if (J-th deriv. at X) <= Y.
                                    !!        = 1,      if (J-th deriv. at X) >= Y.
                                    !!        = 2,      if (J-th deriv. at X) == Y.
                                    !!        = 3,      if (J-th deriv. at X) ==
                                    !!                     (J-th 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 1-4 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 Q-R orthogonal
                                  !! decomposition of the matrix resulting from
                                  !! B-spline 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=NBKPT-NORD, 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 I-th diagonal of R, which has N-I+1
                                  !! elements, is copied to W(*) starting at
                                  !!
                                  !! `W((I-1)*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=NBKPT-NORD
                                     !! parameters are the B-spline 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 B-spline
                                     !! 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=(NBKPT-NORD+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
                    !! B-spline 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 = NBKPT-NORD+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 hard-to-find
                     !! 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,...,I2-1      G(*,*)
   !     I2,...,I3-1      XTEMP(*)
   !     I3,...,I4-1      PTEMP(*)
   !     I4,...,I5-1      BKPT(*) (LOCAL TO [[DFCMN]])
   !     I5,...,I6-1      BF(*,*)
   !     I6,...,I7-1      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*(nbkpt-nord+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