dnedd Subroutine

subroutine dnedd( &
   t, y, ydot, neq, res, jac, dum1, &
   h, wt, jstart, idid, rpar, ipar, phi, gama, dum2, delta, e, &
   rwm, iwm, cj, cjold, cjlast, s, uround, dum3, dum4, dum5, &
   epscon, jcalc, dum6, kp1, nonneg, ntype, iernls)
!! 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(jac_t) :: jac
      !! User-defined Jacobian routine.
   procedure(psol_t) :: dum1
      !! Dummy argument.
   real(rk), intent(in) :: h
      !! Step size.
   real(rk), intent(inout) :: wt(neq) ! @todo: confusing notation: wt ?= ewt ?= whgt
      !! Weights for error control.
   integer, intent(in) :: jstart ! @todo: ? convert to logical
      !! Flag indicating whether this is the first call to this routine.
      !! If `jstart = 0`, then this is the first call, otherwise it is not.
   integer, intent(out) :: idid
      !! Completion flag.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(in) :: phi(neq, kp1)
      !! Array of divided differences used by the Newton iteration.
   real(rk), intent(in) :: gama(kp1)
      !! Array used to predict `y` and `ydot`.
   real(rk), intent(in) :: dum2(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(out) :: cjold
      !! Value of `cj` as of the last call to [[dmatd]]. Accounts for changes in `cj`
      !! needed to decide whether to call [[dmatd]].
   real(rk), intent(in) :: cjlast
      !! Previous value of `cj`.
   real(rk), intent(out) :: s
      !! Convergence factor for the Newton iteration. See [[dnsd]] and [[dnsk]] for
      !! details. On the first Newton iteration with an updated preconditioner, `s = 100`.
      !! The value of `s` is preserved from call to call so that the rate estimate from
      !! a previous step can be applied to the current step.
   real(rk), intent(in) :: uround ! @todo: ? remove
      !! Unit roundoff.
   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.
   integer, intent(out) :: jcalc
      !! Flag indicating whether the Jacobian matrix needs to be updated.
      !! `-1`: call [[dmatd]] to update the Jacobian matrix.
      !!  `0`: Jacobian matrix is up-to-date.
      !!  `1`: Jacobian matrix is outdated, but [[dmatd]] will not be called unless `jcalc` is set to `-1`.
   integer, intent(in) :: dum6
      !! Dummy argument.
   integer, intent(in) :: kp1
      !! The current order + 1. Updated across calls.
   integer, intent(in) :: nonneg
      !! Flag indicating whether nonnegativity constraints are imposed on the solution.
      !! `0`: no nonnegativity constraints.
      !! `1`: nonnegativity constraints are imposed.
   integer, intent(in) :: ntype
      !! Type of the nonlinear solver.
      !! `0`: modified Newton with direct solver.
      !! `1`: modified Newton with iterative linear solver.
      !! `2`: modified Newton with user-supplied linear solver.
   integer, intent(out) :: iernls
      !! Error flag for the nonlinear solver.
      !!  `0`: nonlinear solver converged.
      !!  `1`: recoverable error inside non-linear solver.
      !! `-1`: unrecoverable error inside non-linear solver.

   integer, parameter :: muldel = 1, maxit = 4
   real(rk), parameter :: xrate = 0.25_rk
   integer :: idum, ierj, iernew, iertyp, ires, j
   real(rk) :: delnorm, pnorm, temp1, temp2, tolnew
   real(rk), external :: ddwnrm ! @todo: remove this once inside module

   ! Verify that this is the correct subroutine.
   iertyp = 0
   if (ntype .ne. 0) then
      iertyp = 1
      goto 380
   end if

   ! If this is the first step, perform initializations.
   if (jstart .eq. 0) then
      cjold = cj
      jcalc = -1
   end if

   ! Perform all other initializations.
   iernls = 0

   ! Decide whether new Jacobian is needed.
   temp1 = (one - xrate)/(one + xrate)
   temp2 = 1/temp1
   if ((cj/cjold .lt. temp1) .or. (cj/cjold .gt. temp2)) jcalc = -1
   if (cj .ne. cjlast) s = 1e2_rk

   ! Entry point for updating the Jacobian with current stepsize.
300 continue

   ! Initialize all error flags to zero.
   ierj = 0
   ires = 0
   iernew = 0

   ! Predict the solution and derivative and compute the tolerance for the Newton iteration.
   y = phi(:, 1)
   ydot = zero
   do j = 2, kp1
      y = y + phi(:, j)
      ydot = ydot + gama(j)*phi(:, j)
   end do

   pnorm = ddwnrm(neq, y, wt, rpar, ipar)
   tolnew = 100*uround*pnorm

   ! Call RES to initialize DELTA.
   call res(t, y, ydot, cj, delta, ires, rpar, ipar)
   iwm(iwloc_nre) = iwm(iwloc_nre) + 1
   if (ires .lt. 0) goto 380

   ! If indicated, reevaluate the iteration matrix
   ! J = dG/dY + CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
   ! Set JCALC to 0 as an indicator that this has been done.
   if (jcalc .eq. -1) then
      iwm(iwloc_nje) = iwm(iwloc_nje) + 1
      jcalc = 0
      call dmatd(neq, t, y, ydot, delta, cj, h, ierj, wt, e, rwm, iwm, res, ires, uround, jac, rpar, ipar)
      cjold = cj
      s = 1e2_rk
      if (ires .lt. 0) goto 380
      if (ierj .ne. 0) goto 380
   end if

   ! Call the nonlinear Newton solver.
   temp1 = 2/(one + cj/cjold)
   call dnsd(t, y, ydot, neq, res, dum1, wt, rpar, ipar, dum2, &
             delta, e, rwm, iwm, cj, dum4, dum5, dum3, epscon, s, temp1, &
             tolnew, muldel, maxit, ires, idum, iernew)

   if (iernew .gt. 0 .and. jcalc .ne. 0) then
      ! The Newton iteration had a recoverable failure with an old
      ! iteration matrix. Retry the step with a new iteration matrix.
      jcalc = -1
      goto 300
   end if

   if (iernew .ne. 0) goto 380

   ! The Newton iteration has converged. If nonnegativity of
   ! solution is required, set the solution nonnegative, if the
   ! perturbation to do it is small enough. If the change is too
   ! large, then consider the corrector iteration to have failed.
   if (nonneg .eq. 0) goto 390
   delta = min(y, zero)
   delnorm = ddwnrm(neq, delta, wt, rpar, ipar)
   if (delnorm .gt. epscon) goto 380
   e = e - delta

   ! Exits from nonlinear solver.
   ! No convergence with current iteration matrix, or singular iteration matrix.
   ! Compute IERNLS and IDID accordingly.
380 continue
   if (ires .le. -2 .or. iertyp .ne. 0) then
      iernls = -1
      if (ires .le. -2) idid = -11
      if (iertyp .ne. 0) idid = -15
   else
      iernls = 1
      if (ires .lt. 0) idid = -10
      if (ierj .ne. 0) idid = -8
   end if

390 continue
   jcalc = 1

end subroutine dnedd