dnsk Subroutine

subroutine dnsk( &
   t, y, ydot, neq, res, psol, wt, rpar, ipar, &
   savr, delta, e, rwm, iwm, cj, sqrtn, rsqrtn, epslin, epscon, &
   s, confac, tolnew, muldel, maxit, ires, iersl, iernew)
!! This routines 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.

   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) :: psol
      !! User-defined preconditioner routine.
   real(rk), intent(inout) :: wt(neq) ! @todo: confusing notation: wt ?= ewt ?= whgt
      !! Weights for error control.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(out) :: savr(neq)
      !! Saved residual vector.
   real(rk), intent(inout) :: delta(neq)
      !! Real workspace.
   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) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   real(rk), intent(in) :: epslin
      !! Tolerance for linear system solver.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   real(rk), intent(out) :: s
      !! Convergence factor for the Newton iteration. `s=rate/(1 - rate)`, where
      !! `rate` is the estimated rate of convergence of the Newton iteration.
      !! The closer `rate` is to 0, the faster the Newton iteration is converging.
      !! the closer `rate` is to 1, the slower the Newton iteration is converging.
   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`.
   integer, intent(out) :: iersl
      !! Error flag from the linear system solver. See [[dslvk]] for details.
      !! If `iersl < 0`, then `iernew = -1`.
      !! If `iersl == 1`, then `iernew = 1`.
   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, rhok
   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 == 1) then
         delta = delta*confac
      end if

      ! Save residual in SAVR.
      savr = delta

      ! Compute a new iterate. Store the correction in DELTA.
      call dslvk(neq, y, t, ydot, savr, delta, wt, rwm, iwm, &
                 res, ires, psol, iersl, cj, epslin, sqrtn, rsqrtn, rhok, &
                 rpar, ipar)
      if ((ires /= 0) .or. (iersl /= 0)) exit

      ! 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 == 0) then
         oldnorm = deltanorm
         if (deltanorm <= tolnew) then
            converged = .true.
            exit
         end if
      else
         rate = (deltanorm/oldnorm)**(one/m)
         if (rate > 0.9_rk) exit
         s = rate/(one - rate)
      end if

      if (s*deltanorm <= 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 >= 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 < 0) exit

   end do

   if (converged) then
      iernew = 0
   else
      if ((ires <= -2) .or. (iersl < 0)) then
         iernew = -1
      else
         iernew = 1
      end if
   end if

end subroutine dnsk