dnsd Subroutine

subroutine dnsd( &
   t, y, ydot, neq, res, dum1, wt, rpar, ipar, &
   dum2, delta, e, rwm, iwm, cj, dum3, dum4, dum5, epscon, &
   s, confac, tolnew, muldel, maxit, ires, dum6, iernew)
!! This routine solves a nonlinear system of algebraic equations of the form:
!!
!!  $$ G(t, y, \dot{y}) = 0 $$
!!
!! for the unknown \(y\). The method used is a modified Newton scheme.
!!
!! All arguments with "dum" in their names are dummy arguments which are not used in
!! this routine.

   use daskr_kinds, only: rk, zero, one
   use daskr
   implicit none

   real(rk), intent(in) :: t
      !! Independent variable.
   real(rk), intent(inout) :: y(neq)
      !! Solution vector.
   real(rk), intent(inout) :: ydot(neq)
      !! Derivative of solution vector.
   integer, intent(in) :: neq
      !! Problem size.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(psol_t) :: dum1
      !! Dummy argument.
   real(rk), intent(inout) :: wt(neq)
      !! Weights for error criterion.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(in) :: dum2(neq)
      !! Dummy argument.
   real(rk), intent(inout) :: delta(neq)
      !! Real work array.
   real(rk), intent(out) :: e(neq)
      !! Error accumulation vector.
   real(rk), intent(inout) :: rwm(*)
      !! Real workspace for the linear system solver.
   integer, intent(inout) :: iwm(*)
      !! Integer workspace for the linear system solver.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: dum3
      !! Dummy argument.
   real(rk), intent(in) :: dum4
      !! Dummy argument.
   real(rk), intent(in) :: dum5
      !! Dummy argument.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   real(rk), intent(inout) :: s
      !! Convergence factor for the Newton iteration. `s=rate/(1 - rate)`, where
      !! `rate` is the estimated rate of convergence of the Newton iteration.
   real(rk), intent(in) :: confac
      !! Residual scale factor to improve convergence.
   real(rk), intent(in) :: tolnew
      !! Tolerance on the norm of the Newton correction in the alternative Newton
      !! convergence test.
   integer, intent(in) :: muldel ! @todo: convert to logical
      !! Flag indicating whether or not to multiply `delta` by `confac`.
   integer, intent(in) :: maxit
      !! Maximum number of Newton iterations allowed.
   integer, intent(out) :: ires
      !! Error flag from `res`.
      !! If `ires == -1`, then `iernew = 1`.
      !! If `ires < -1`, then `iernew = -1`.
   integer, intent(in) :: dum6
      !! Dummy argument.
   integer, intent(out) :: iernew
      !! Error flag for the Newton iteration.
      !!  `0`: the Newton iteration converged.
      !!  `1`: recoverable error inside Newton iteration.
      !! `-1`: unrecoverable error inside Newton iteration.

   integer :: m
   real(rk) :: deltanorm, oldnorm, rate
   real(rk), external :: ddwnrm ! @todo: remove this once inside module
   logical :: converged

   ! Initialize Newton counter M and accumulation vector E.
   m = 0
   e = zero

   ! Corrector loop.
   converged = .false.
   do

      iwm(iwloc_nni) = iwm(iwloc_nni) + 1

      ! If necessary, multiply residual by convergence factor.
      if (muldel .eq. 1) then
         delta = delta*confac
      end if

      ! Compute a new iterate (back-substitution).
      ! Store the correction in DELTA.
      call dslvd(neq, delta, rwm, iwm)

      ! Update Y, E, and YDOT.
      y = y - delta
      e = e - delta
      ydot = ydot - cj*delta

      ! Test for convergence of the iteration.
      deltanorm = ddwnrm(neq, delta, wt, rpar, ipar)
      if (m .eq. 0) then
         oldnorm = deltanorm
         if (deltanorm .le. tolnew) then
            converged = .true.
            exit
         end if
      else
         rate = (deltanorm/oldnorm)**(one/m)
         if (rate .gt. 0.9_rk) exit
         s = rate/(1 - rate)
      end if
      if (s*deltanorm .le. epscon) then
         converged = .true.
         exit
      end if

      ! The corrector has not yet converged.
      ! Update M and test whether the maximum number of iterations have been tried.
      m = m + 1
      if (m .ge. maxit) exit

      ! Evaluate the residual, and go back to do another iteration.
      ires = 0
      call res(t, y, ydot, cj, delta, ires, rpar, ipar)
      iwm(iwloc_nre) = iwm(iwloc_nre) + 1
      if (ires .lt. 0) exit

   end do

   if (converged) then
      iernew = 0
   else
      if (ires .le. -2) then
         iernew = -1
      else
         iernew = 1
      end if
   end if

end subroutine dnsd