defcmn Subroutine

   subroutine defcmn(Ndata, Xdata, Ydata, Sddata, Nord, Nbkpt, Bkptin, &
                     Mdein, Mdeout, Coeff, Bf, Xtemp, Ptemp, Bkpt, g, Mdg, w, &
                     Mdw, Lw)

      integer :: Lw, Mdein, Mdeout, Mdg, Mdw, Nbkpt, Ndata, Nord
      real(wp) :: Bf(Nord, *), Bkpt(*), Bkptin(*), Coeff(*), &
                  g(Mdg, *), Ptemp(*), Sddata(*), w(Mdw, *), &
                  Xdata(*), Xtemp(*), Ydata(*)

      real(wp) :: rnorm, xmax, xmin, xval
      integer :: i, idata, ileft, intseq, ip, ir, irow, l, mt, n, &
                 nb, nordm1, nordp1, np1
      character(len=8) :: xern1, xern2
      integer :: dfspvn_j
      real(wp), dimension(20) :: dfspvn_deltam, dfspvn_deltap

      ! Initialize variables and analyze input.

      n = Nbkpt - Nord
      np1 = n + 1

      ! Initially set all output coefficients to zero.

      call dcopy(n, [0.0_wp], 0, Coeff, 1)

      Mdeout = -1
      if (Nord < 1 .or. Nord > 20) then
         write (*, *) 'IN DEFC, THE ORDER OF THE B-SPLINE MUST BE 1 THRU 20.'
         return
      end if

      if (Nbkpt < 2*Nord) then
         write (*, *) 'IN DEFC, THE NUMBER OF KNOTS MUST BE AT LEAST TWICE THE B-SPLINE ORDER.'
         return
      end if

      if (Ndata < 0) then
         write (*, *) 'IN DEFC, THE NUMBER OF DATA POINTS MUST BE NONNEGATIVE.'
         return
      end if

      nb = (Nbkpt - Nord + 3)*(Nord + 1) + (Nbkpt + 1)*(Nord + 1) &
           + 2*max(Nbkpt, Ndata) + Nbkpt + Nord**2
      if (Lw < nb) then
         write (xern1, '(I8)') nb
         write (xern2, '(I8)') Lw
         write (*, *) 'IN DEFC, INSUFFICIENT STORAGE FOR W(*).  CHECK FORMULA '// &
            'THAT READS LW>= ... .  NEED = '//xern1// &
            ' GIVEN = '//xern2
         Mdeout = -1
         return
      end if

      if (Mdein /= 1 .and. Mdein /= 2) then
         write (*, *) 'IN DEFC, INPUT VALUE OF MDEIN MUST BE 1-2.'
         return
      end if

      ! Sort the breakpoints.

      call dcopy(Nbkpt, Bkptin, 1, Bkpt, 1)
      call dsort(Nbkpt, 1, Bkpt)

      ! Save interval containing knots.

      xmin = Bkpt(Nord)
      xmax = Bkpt(np1)
      nordm1 = Nord - 1
      nordp1 = Nord + 1

      ! 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
      ! JW : really Ptemp should be an integer array.
      !      it is real because they are stuffing it in
      !      a real work array and also using dsort on it.

      if (Ndata > 0) then
         call dsort(Ndata, 2, Xtemp, Ptemp)
         xmin = min(xmin, Xtemp(1))
         xmax = max(xmax, Xtemp(Ndata))
      end if

      ! Fix breakpoint array if needed. This should only involve very
      ! minor differences with the input array of breakpoints.

      do i = 1, Nord
         Bkpt(i) = min(Bkpt(i), xmin)
      end do

      do i = np1, Nbkpt
         Bkpt(i) = max(Bkpt(i), xmax)
      end do

      ! set up variables for dfspvn
      dfspvn_j = 1
      dfspvn_deltam = 0.0_wp
      dfspvn_deltap = 0.0_wp

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

      mt = 0
      ip = 1
      ir = 1
      ileft = Nord
      intseq = 1
      do idata = 1, Ndata

         ! Sorted indices are in PTEMP(*).

         l = int(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<=N.

            do ileft = ileft, n
               if (xval < Bkpt(ileft + 1)) exit
               if (Mdein == 2) then
                  !  Data is being sequentially accumulated.
                  !  Transfer previously accumulated rows from W(*,*) to
                  !  G(*,*) and process them.
                  call dcopy(nordp1, w(intseq, 1), Mdw, g(ir, 1), Mdg)
                  call dbndac(g, Mdg, Nord, ip, ir, 1, intseq)
                  intseq = intseq + 1
               end if
            end do
         end if

         ! 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
         end if
      end do

      ! Process last block of equations.

      call dbndac(g, Mdg, Nord, ip, ir, mt, ileft - nordm1)

      ! Finish processing any previously accumulated rows from W(*,*)
      ! to G(*,*).

      if (Mdein == 2) then
         do i = intseq, np1
            call dcopy(nordp1, w(i, 1), Mdw, g(ir, 1), Mdg)
            call dbndac(g, Mdg, Nord, ip, ir, 1, min(n, i))
         end do
      end if

      ! 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)

      ! Transfer accumulated rows from G(*,*) to W(*,*) for
      ! possible later sequential accumulation.

      do i = 1, np1
         call dcopy(nordp1, g(i, 1), Mdg, w(i, 1), Mdw)
      end do

     ! Solve for coefficients when possible.

      do i = 1, n
         if (g(i, 1) == 0.0_wp) then
            Mdeout = 2
            return
         end if
      end do

      ! All the diagonal terms in the accumulated triangular
      ! matrix are nonzero.  The solution can be computed but
      ! it may be unsuitable for further use due to poor
      ! conditioning or the lack of constraints.  No checking
      ! for either of these is done here.

      call dbndsl(1, g, Mdg, Nord, ip, ir, Coeff, n, rnorm)
      Mdeout = 1

   end subroutine defcmn