dfcmn Subroutine

   subroutine dfcmn (ndata, xdata, ydata, sddata, nord, nbkpt, &
                     bkptin, nconst, xconst, yconst, nderiv, mode, coeff, bf, xtemp, &
                     ptemp, bkpt, g, mdg, w, mdw, work, iwork)

   integer :: iwork(*), mdg, mdw, mode, nbkpt, nconst, ndata, nderiv(*), &
              nord
   real(wp) :: bf(nord,*), bkpt(*), bkptin(*), coeff(*), &
               g(mdg,*), ptemp(*), sddata(*), w(mdw,*), work(*), &
               xconst(*), xdata(*), xtemp(*), yconst(*), ydata(*)

   real(wp) :: prgopt(10), rnorm, rnorme, rnorml, xmax, &
               xmin, xval, yval
   integer :: i, idata, ideriv, ileft, intrvl, intw1, ip, ir, irow, &
              itype, iw1, iw2, l, lw, mt, n, nb, neqcon, nincon, nordm1, &
              nordp1, np1
   logical :: band, new, var
   character(len=8) :: xern1
   integer :: dfspvn_j
   real(wp), dimension(20) :: dfspvn_deltam, dfspvn_deltap

   ! Analyze input.

   if (nord<1 .or. nord>20) then
      write(*,*) 'IN DFC, THE ORDER OF THE B-SPLINE MUST BE 1 THRU 20.'
      mode = -1
      return
   elseif (nbkpt<2*nord) then
      write(*,*) 'IN DFC, THE NUMBER OF KNOTS MUST BE AT LEAST TWICE ' // &
                 'THE B-SPLINE ORDER.'
      mode = -1
      return
   endif
   if (ndata<0) then
      write(*,*) 'IN DFC, THE NUMBER OF DATA POINTS MUST BE NONNEGATIVE.'
      mode = -1
      return
   endif

   ! Amount of storage allocated for W(*), IW(*).

   iw1 = iwork(1)
   iw2 = iwork(2)
   nb = (nbkpt-nord+3)*(nord+1) + 2*max(ndata,nbkpt) + nbkpt + &
        nord**2

   ! See if sufficient storage has been allocated.

   if (iw1<nb) then
      write (xern1, '(I8)') nb
      write(*,*) 'IN DFC, INSUFFICIENT STORAGE FOR W(*).  CHECK NB = ' // xern1
      mode = -1
      return
   endif

   select case (mode)
   case (1)
      band = .true.
      var = .false.
      new = .true.
   case (2)
      band = .false.
      var = .true.
      new = .true.
   case (3)
      band = .true.
      var = .false.
      new = .false.
   case (4)
      band = .false.
      var = .true.
      new = .false.
   case default
      write(*,*) 'IN DFC, INPUT VALUE OF MODE MUST BE 1-4.'
      mode = -1
      return
   end select
   mode = 0

   ! Sort the breakpoints.

   call dcopy (nbkpt, bkptin, 1, bkpt, 1)
   call dsort (nbkpt, 1, bkpt)

   ! Initialize variables.

   neqcon = 0
   nincon = 0
   do i = 1,nconst
      l = nderiv(i)
      itype = mod(l,4)
      if (itype<2) then
         nincon = nincon + 1
      else
         neqcon = neqcon + 1
      endif
   end do
   ! set up variables for dfspvn
   dfspvn_j = 1
   dfspvn_deltam = 0.0_wp
   dfspvn_deltap = 0.0_wp

   ! Compute the number of variables.

   n = nbkpt - nord
   np1 = n + 1
   lw = nb + (np1+nconst)*np1 + 2*(neqcon+np1) + (nincon+np1) + &
        (nincon+2)*(np1+6)
   intw1 = nincon + 2*np1

   ! Save interval containing knots.

   xmin = bkpt(nord)
   xmax = bkpt(np1)

   ! Find the smallest referenced independent variable value in any
   ! constraint.

   do i = 1,nconst
      xmin = min(xmin,xconst(i))
      xmax = max(xmax,xconst(i))
   end do
   nordm1 = nord - 1
   nordp1 = nord + 1

   ! Define the option vector PRGOPT(1-10) for use in [[DLSEI]].

   prgopt(1) = 4

   ! Set the covariance matrix computation flag.

   prgopt(2) = 1
   if (var) then
      prgopt(3) = 1
   else
      prgopt(3) = 0
   endif

   ! Increase the rank determination tolerances for both equality
   ! constraint equations and least squares equations.

   prgopt(4) = 7
   prgopt(5) = 4
   prgopt(6) = 1.0e-4_wp
   prgopt(7) = 10
   prgopt(8) = 5
   prgopt(9) = 1.0e-4_wp
   prgopt(10) = 1

   ! Turn off work array length checking in [[DLSEI]].

   iwork(1) = 0
   iwork(2) = 0

   ! Initialize variables and analyze input.

   if (new) then

      ! To process least squares equations sort data and an array of
      ! pointers.

      call dcopy (ndata, xdata, 1, xtemp, 1)
      do i = 1,ndata
         ptemp(i) = i
      end do

      if (ndata>0) then
         call dsort (ndata, 2, xtemp, ptemp)
         xmin = min(xmin,xtemp(1))
         xmax = max(xmax,xtemp(ndata))
      endif

      ! Fix breakpoint array if needed.

      do i = 1,nord
         bkpt(i) = min(bkpt(i),xmin)
      end do

      do i = np1,nbkpt
         bkpt(i) = max(bkpt(i),xmax)
      end do

      ! Initialize parameters of banded matrix processor, DBNDAC( ).

      mt = 0
      ip = 1
      ir = 1
      ileft = nord
      do idata = 1,ndata

         ! Sorted indices are in PTEMP(*).

         l = ptemp(idata)
         xval = xdata(l)

         ! When interval changes, process equations in the last block.

         if (xval>=bkpt(ileft+1)) then
            call dbndac (g, mdg, nord, ip, ir, mt, ileft-nordm1)
            mt = 0

            ! Move pointer up to have BKPT(ILEFT)<=XVAL,
            ! ILEFT<NP1.

            do
               if (xval>=bkpt(ileft+1) .and. ileft<n) then
                  ileft = ileft + 1
               else
                  exit
               endif
            end do

         endif

         ! Obtain B-spline function value.

         call dfspvn (bkpt, nord, 1, xval, ileft, bf, &
                      dfspvn_j, dfspvn_deltam, dfspvn_deltap)

         ! Move row into place.

         irow = ir + mt
         mt = mt + 1
         call dcopy (nord, bf, 1, g(irow,1), mdg)
         g(irow,nordp1) = ydata(l)

         ! Scale data if uncertainty is nonzero.

         if (sddata(l)/=0.0_wp) call dscal (nordp1, 1.0_wp/sddata(l), &
                                            g(irow,1), mdg)

         ! When staging work area is exhausted, process rows.

         if (irow==mdg-1) then
            call dbndac (g, mdg, nord, ip, ir, mt, ileft-nordm1)
            mt = 0
         endif
      end do

      ! Process last block of equations.

      call dbndac (g, mdg, nord, ip, ir, mt, ileft-nordm1)

      ! Last call to adjust block positioning.

      call dcopy (nordp1, [0.0_wp], 0, g(ir,1), mdg)
      call dbndac (g, mdg, nord, ip, ir, 1, np1)
   endif

   band = band .and. nconst==0
   do i = 1,n
      band = band .and. g(i,1)/=0.0_wp
   end do

   ! Process banded least squares equations.

   if (band) then
      call dbndsl (1, g, mdg, nord, ip, ir, coeff, n, rnorm)
      return
   endif

   ! Check further for sufficient storage in working arrays.

   if (iw1<lw) then
      write (xern1, '(I8)') lw
      write(*,*) 'IN DFC, INSUFFICIENT STORAGE FOR W(*).  CHECK LW = ' // xern1
      mode = -1
      return
   endif

   if (iw2<intw1) then
      write (xern1, '(I8)') intw1
      write(*,*) 'IN DFC, INSUFFICIENT STORAGE FOR IW(*).  CHECK IW1 = ' // xern1
      mode = -1
      return
   endif

   ! Write equality constraints.
   ! Analyze constraint indicators for an equality constraint.

   neqcon = 0
   do idata = 1,nconst
      l = nderiv(idata)
      itype = mod(l,4)
      if (itype>1) then
         ideriv = l/4
         neqcon = neqcon + 1
         ileft = nord
         xval = xconst(idata)

         do
            if (xval<bkpt(ileft+1) .or. ileft>=n) exit
            ileft = ileft + 1
         end do

         call dfspvd (bkpt, nord, xval, ileft, bf, ideriv+1)
         call dcopy (np1, [0.0_wp], 0, w(neqcon,1), mdw)
         call dcopy (nord, bf(1,ideriv+1), 1, w(neqcon,ileft-nordm1), &
                     mdw)

         if (itype==2) then
            w(neqcon,np1) = yconst(idata)
         else
            ileft = nord
            yval = yconst(idata)

            do
               if (yval<bkpt(ileft+1) .or. ileft>=n) exit
               ileft = ileft + 1
            end do

            call dfspvd (bkpt, nord, yval, ileft, bf, ideriv+1)
            call daxpy (nord, -1.0_wp, bf(1, ideriv+1), 1, &
                        w(neqcon, ileft-nordm1), mdw)
         endif
      endif
   end do

   ! Transfer least squares data.

   do i = 1,np1
      irow = i + neqcon
      call dcopy (n, [0.0_wp], 0, w(irow,1), mdw)
      call dcopy (min(np1-i, nord), g(i,1), mdg, w(irow,i), mdw)
      w(irow,np1) = g(i,nordp1)
   end do

   ! Write inequality constraints.
   ! Analyze constraint indicators for inequality constraints.

   nincon = 0
   do idata = 1,nconst
      l = nderiv(idata)
      itype = mod(l,4)
      if (itype<2) then
         ideriv = l/4
         nincon = nincon + 1
         ileft = nord
         xval = xconst(idata)

         do
            if (xval<bkpt(ileft+1) .or. ileft>=n) exit
            ileft = ileft + 1
         end do

         call dfspvd (bkpt, nord, xval, ileft, bf, ideriv+1)
         irow = neqcon + np1 + nincon
         call dcopy (n, [0.0_wp], 0, w(irow,1), mdw)
         intrvl = ileft - nordm1
         call dcopy (nord, bf(1, ideriv+1), 1, w(irow, intrvl), mdw)

         if (itype==1) then
            w(irow,np1) = yconst(idata)
         else
            w(irow,np1) = -yconst(idata)
            call dscal (nord, -1.0_wp, w(irow, intrvl), mdw)
         endif
      endif
   end do

   ! Solve constrained least squares equations.

   call dlsei(w, mdw, neqcon, np1, nincon, n, prgopt, coeff, rnorme, &
              rnorml, mode, work, iwork)

   end subroutine dfcmn