daskr_new.f90 Source File

!----------------------------------------------------------------------------------------------
! Adapted from original Fortran code in `original/solver/ddaskr.f`
!----------------------------------------------------------------------------------------------

module daskr

   use daskr_kinds, only: rk
   implicit none

   abstract interface
      subroutine res_t(t, y, ydot, cj, delta, ires, rpar, ipar)
         import :: rk
         real(rk), intent(in) :: t
         real(rk), intent(in) :: y(*)
         real(rk), intent(in) :: ydot(*)
         real(rk), intent(in) :: cj
         real(rk), intent(out) :: delta(*)
         integer, intent(inout) :: ires
         real(rk), intent(inout) :: rpar(*)
         integer, intent(inout) :: ipar(*)
      end subroutine res_t

      subroutine jac_t(t, y, ydot, pd, cj, rpar, ipar)
         import :: rk
         procedure(res_t) :: res
         real(rk), intent(in) :: t
         real(rk), intent(in) :: y(*)
         real(rk), intent(in) :: ydot(*)
         real(rk), intent(out) :: pd(*)
         real(rk), intent(in) :: cj
         real(rk), intent(inout) :: rpar(*)
         integer, intent(inout) :: ipar(*)
      end subroutine jac_t

      subroutine jack_t(res, ires, neq, t, y, ydot, rewt, savr, wk, h, cj, rwp, iwp, ierr, rpar, ipar)
         import :: rk
         procedure(res_t) :: res
         integer, intent(out) :: ires
         integer, intent(in) :: neq
         real(rk), intent(in) :: t
         real(rk), intent(inout) :: y(*)
         real(rk), intent(inout) :: ydot(*)
         real(rk), intent(in) :: rewt(*)
         real(rk), intent(inout) :: savr(*)
         real(rk), intent(inout) :: wk(*) !(neq)
         real(rk), intent(in) :: h
         real(rk), intent(in) :: cj
         real(rk), intent(inout) :: rwp(*)
         integer, intent(inout) :: iwp(*)
         integer, intent(inout) :: ierr
         real(rk), intent(inout) :: rpar(*)
         integer, intent(inout) :: ipar(*)
      end subroutine jack_t

      subroutine psol_t(neq, t, y, ydot, savr, wk, cj, wght, rwp, iwp, b, epslin, ierr, rpar, ipar)
         import :: rk
         integer, intent(in) :: neq
         real(rk), intent(in) :: t
         real(rk), intent(in) :: y(*)
         real(rk), intent(in) :: ydot(*)
         real(rk), intent(in) :: savr(*)
         real(rk), intent(inout) :: wk(*)
         real(rk), intent(in) :: cj
         real(rk), intent(in) :: wght(*)
         real(rk), intent(inout) :: rwp(*)
         integer, intent(inout) :: iwp(*)
         real(rk), intent(inout) :: b(*)
         real(rk), intent(in) :: epslin
         integer, intent(inout) :: ierr
         real(rk), intent(inout) :: rpar(*)
         integer, intent(inout) :: ipar(*)
      end subroutine psol_t
   end interface

   ! `rwork` indices
   integer, parameter :: rwloc_cj = 1
      !! `cj`: scalar used in forming the system Jacobian
   integer, parameter :: rwloc_cjold = 2
      !! `cjold`: scalar used in forming the system Jacobian on last step
   integer, parameter :: rwloc_h = 3
      !! `h`: step size for next step
   integer, parameter :: rwloc_tn = 4
      !! `tn`: current independent variable value
   integer, parameter :: rwloc_hold = 7
      !! `hold`: step size used on last step

   ! `iwork` indices
   integer, parameter :: iwloc_ml = 1
      !! `ml`: lower bandwidth
   integer, parameter :: iwloc_mu = 2
      !! `mu`: upper bandwidth
   integer, parameter :: iwloc_mxord = 3
      !! `mxord`: method order
   integer, parameter :: iwloc_mtype = 4
      !! `mtype`: method type
   integer, parameter :: iwloc_k = 7
      !! `k`: order of method for next step
   integer, parameter :: iwloc_kold = 8
      !! `kold`: order used on last step
   integer, parameter :: iwloc_nst = 11
      !! `nst`: number of steps taken
   integer, parameter :: iwloc_nre = 12
      !! `nre`: number of `res` calls
   integer, parameter :: iwloc_nje = 13
      !! `nje`: number of `jac` calls
   integer, parameter :: iwloc_netf = 14
      !! `netf`: number of error test failures
   integer, parameter :: iwloc_ncfn = 15
      !! `ncfn`: number of nonlinear convergence failures
   integer, parameter :: iwloc_ncfl = 16
      !! `ncfl`: number of linear iteration convergence failures
   integer, parameter :: iwloc_leniw = 17
      !! `leniw`: actual `iwork` length used
   integer, parameter :: iwloc_lenrw = 18
      !! `lenrw`: actual `rwork` length used
   integer, parameter :: iwloc_nni = 19
      !! `nni`: number of nonlinear iterations
   integer, parameter :: iwloc_nli = 20
      !! `nli`: number of linear (Krylov) iterations
   integer, parameter :: iwloc_nps = 21
      !! `nps`: number of `psol` calls
   integer, parameter :: iwloc_npd = 22
      !! `npd`: length of partial derivatives vector returned by `jac`
   integer, parameter :: iwloc_miter = 23
      !! `miter`: ??
   integer, parameter :: iwloc_maxl = 24
      !! `maxl`: ??
   integer, parameter :: iwloc_kmp = 25
      !! `kmp`: ??
   integer, parameter :: iwloc_nrmax = 26
      !! `nrmax`: ??
   integer, parameter :: iwloc_lrwp = 29
      !! `lrwp`: location of start of `rwp` in `rwork`.
   integer, parameter :: iwloc_liwp = 30
      !! `liwp`: location of start of `iwp` in `iwork`.
   integer, parameter :: iwloc_kprint = 31
      !! `kprint`: flag to turn on print
   integer, parameter :: iwloc_mxnit = 32
      !! `mxnit`: maximum number of Newton iterations
   integer, parameter :: iwloc_mxnj = 33
      !! `mxnj`: maximum number of Jacobian evaluations ?
   integer, parameter :: iwloc_lsoff = 35
      !! `lsoff`: flag to turn off linesearch
   integer, parameter :: iwloc_nrtfn = 36
      !! `nrtfn`: number of `rt` calls

end module daskr

subroutine ddstp( &
   t, y, ydot, neq, res, jac, psol, h, wt, vt, &
   jstart, idid, rpar, ipar, phi, savr, delta, e, rwm, iwm, &
   alpha, beta, gama, psi, sigma, cj, cjold, hold, s, hmin, uround, &
   epli, sqrtn, rsqrtn, epscon, iphase, jcalc, jflg, k, kold, ns, nonneg, ntype, nls)

!! This routine solves a system of differential/algebraic equations of the form:
!!
!! $$ G(t,y,\dot{y}) = 0 $$
!!
!! for one step (normally from \(t\) to \(t+h\)). The methods used are modified divided
!! difference, fixed leading coefficient forms of backward differentiation formulas.
!! The code adjusts the stepsize and order to control the local error per step.

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

   real(rk), intent(inout) :: t
      !! Independent variable.
   real(rk), intent(inout) :: y(neq)
      !! Solution vector.
   real(rk), intent(inout) :: ydot(neq)
      !! Derivative of solution vector after successful step.
   integer, intent(in) :: neq
      !! Problem size.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(jack_t) :: jac
      !! User-defined Jacobian routine.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   real(rk), intent(inout) :: h
      !! Step size.
   real(rk), intent(inout) :: wt(neq)
      !! Weights for error control.
   real(rk), intent(inout) :: vt(neq)
      !! Masked weights for error control.
   integer, intent(inout) :: 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(inout) :: phi(neq, *)
      !! Array of divided differences used by the Newton iteration.
   real(rk), intent(out) :: savr(neq)
      !! Real work array.
   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(out) :: alpha(*)
      !! ??
   real(rk), intent(out) :: beta(*)
      !! ??
   real(rk), intent(out) :: gama(*)
      !! ??
   real(rk), intent(inout) :: psi(*)
      !! ??
   real(rk), intent(out) :: sigma(*)
      !! ??
   real(rk), intent(inout) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(out) :: cjold
      !! Value of `cj` as of the last call to [[ditmd]]. Accounts for changes in `cj`
      !! needed to decide whether to call [[ditmd]].
   real(rk), intent(out) :: hold
      !! ??
   real(rk), intent(out) :: s
      !! Convergence factor for the Newton iteration.
   real(rk), intent(in) :: hmin
      !! ??
   real(rk), intent(in) :: uround
      !! Unit roundoff.
   real(rk), intent(in) :: epli ! @note: not the same as 'epslin'!
      !! Convergence test constant for the iterative solution of linear systems. See
      !! [[daskr]] for details.
   real(rk), intent(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   integer, intent(inout) :: iphase
      !! ??
   integer, intent(out) :: jcalc
      !! Flag indicating whether the Jacobian matrix needs to be updated.
   integer, intent(in) :: jflg
      !! Flag indicating whether a `jac` routine is supplied by the user.
   integer, intent(inout) :: k
      !! ??
   integer, intent(inout) :: kold
      !! ??
   integer, intent(inout) :: ns
      !! ??
   integer, intent(in) :: nonneg
      !! Flag indicating whether nonnegativity constraints are imposed on the solution.
   integer, intent(in) :: ntype
      !! Type of the nonlinear solver.
   procedure() :: nls
      !! ??

   integer :: i, iernls, j, j1, kdiff, km1, knew, kp1, kp2, ncf, nef, nsp1
   real(rk) :: alpha0, alphas, cjlast, ck, enorm, erk, erkm1, erkm2, erkp1, err, est, hnew, &
               r, temp1, temp2, terk, terkm1, terkm2, terkp1, told
   real(rk), external :: ddwnrm ! @todo: remove this once inside module

   ! BLOCK 1.
   ! Initialize. On the first call, set the order to 1 and initialize other variables.

   ! Initializations for all calls
   told = t
   ncf = 0
   nef = 0

   ! If this is the first step, perform other initializations
   if (jstart == 0) then
      k = 1
      kold = 0
      hold = zero
      psi(1) = h
      cj = 1/h
      iphase = 0
      ns = 0
   end if

   ! BLOCK 2
   ! Compute coefficients of formulas for this step.
   do

      kp1 = k + 1
      kp2 = k + 2
      km1 = k - 1
      if ((h /= hold) .or. (k /= kold)) ns = 0
      ns = min(ns + 1, kold + 2)
      nsp1 = ns + 1

      if (kp1 >= ns) then

         beta(1) = one
         alpha(1) = one
         temp1 = h
         gama(1) = zero
         sigma(1) = one

         do i = 2, kp1
            temp2 = psi(i - 1)
            psi(i - 1) = temp1
            beta(i) = beta(i - 1)*psi(i - 1)/temp2
            temp1 = temp2 + h
            alpha(i) = h/temp1
            sigma(i) = (i - 1)*sigma(i - 1)*alpha(i)
            gama(i) = gama(i - 1) + alpha(i - 1)/h
         end do

         psi(kp1) = temp1

      end if

      ! Compute ALPHAS, ALPHA0
      alphas = zero
      alpha0 = zero
      do i = 1, k
         alphas = alphas - one/i
         alpha0 = alpha0 - alpha(i)
      end do

      ! Compute leading coefficient CJ
      cjlast = cj
      cj = -alphas/h

      ! Compute variable stepsize error coefficient CK
      ck = abs(alpha(kp1) + alphas - alpha0)
      ck = max(ck, alpha(kp1))

      ! Change PHI to PHI STAR
      if (kp1 >= nsp1) then
         do j = nsp1, kp1
            phi(:, j) = beta(j)*phi(:, j)
         end do
      end if

      ! Update time
      t = t + h

      ! Initialize IDID to 1
      idid = 1

      ! BLOCK 3
      ! Call the nonlinear system solver to obtain the solution and derivative.
      call nls(t, y, ydot, neq, &
               res, jac, psol, h, wt, jstart, idid, rpar, ipar, phi, gama, &
               savr, delta, e, rwm, iwm, cj, cjold, cjlast, s, &
               uround, epli, sqrtn, rsqrtn, epscon, jcalc, jflg, kp1, &
               nonneg, ntype, iernls)

      if (iernls /= 0) goto 600

      ! BLOCK 4
      ! Estimate the errors at orders K,K-1,K-2 as if constant stepsize was used.
      ! Estimate the local error at order K and test whether the current step is successful.

      ! Estimate errors at orders K,K-1,K-2
      enorm = ddwnrm(neq, e, vt, rpar, ipar)
      erk = sigma(k + 1)*enorm
      terk = (k + 1)*erk
      est = erk
      knew = k

      if (k == 1) goto 430
      delta = phi(:, kp1) + e
      erkm1 = sigma(k)*ddwnrm(neq, delta, vt, rpar, ipar)
      terkm1 = k*erkm1
      if (k > 2) goto 410
      if (terkm1 <= terk/2) goto 420
      goto 430

410   continue
      delta = delta + phi(:, k)
      erkm2 = sigma(k - 1)*ddwnrm(neq, delta, vt, rpar, ipar)
      terkm2 = (k - 1)*erkm2
      if (max(terkm1, terkm2) > terk) goto 430

      ! Lower the order
420   continue
      knew = k - 1
      est = erkm1

      ! Calculate the local error for the current step to see if the step was successful
430   continue
      err = ck*enorm
      if (err > one) goto 600

      ! BLOCK 5
      ! The step is successful. Determine the best order and stepsize for
      ! the next step. Update the differences for the next step.
      idid = 1
      iwm(iwloc_nst) = iwm(iwloc_nst) + 1
      kdiff = k - kold
      kold = k
      hold = h

      ! Estimate the error at order K+1 unless already decided to lower order,
      ! or already using maximum order,
      ! or stepsize not constant,
      ! or order raised in previous step
      if (knew == km1 .or. k == iwm(iwloc_mxord)) iphase = 1
      if (iphase == 0) goto 545
      if (knew == km1) goto 540
      if (k == iwm(iwloc_mxord)) goto 550
      if (kp1 >= ns .or. kdiff == 1) goto 550

      delta = e - phi(:, kp2)
      erkp1 = ddwnrm(neq, delta, vt, rpar, ipar)/(k + 2)
      terkp1 = (k + 2)*erkp1
      if (k > 1) goto 520
      if (terkp1 >= terk/2) goto 550
      goto 530

520   continue
      if (terkm1 <= min(terk, terkp1)) goto 540
      if (terkp1 >= terk .or. k == iwm(iwloc_mxord)) goto 550

      ! Raise order
530   continue
      k = kp1
      est = erkp1
      goto 550

      ! Lower order
540   continue
      k = km1
      est = erkm1
      goto 550

      ! If IPHASE = 0, increase order by one and multiply stepsize by factor two
545   continue
      k = kp1
      hnew = 2*h
      h = hnew
      goto 575

      ! Determine the appropriate stepsize for the next step.
550   continue
      hnew = h
      temp2 = real(k + 1, rk)
      r = (2*est + 1e-4_rk)**(-1/temp2)
      if (r < two) goto 555
      hnew = 2*h
      goto 560

555   continue
      if (r > one) goto 560
      r = max(0.5_rk, min(0.9_rk, r))
      hnew = h*r

560   continue
      h = hnew

      ! Update differences for next step
575   continue
      if (kold == iwm(iwloc_mxord)) goto 585
      phi(:, kp2) = e

585   continue
      phi(:, kp1) = phi(:, kp1) + e
      do j1 = 2, kp1
         j = kp1 - j1 + 1
         phi(:, j) = phi(:, j) + phi(:, j + 1)
      end do
      jstart = 1
      return

      ! BLOCK 6
      ! The step is unsuccessful. Restore X, PSI, PHI
      ! Determine appropriate stepsize for continuing the integration,
      ! or exit with an error flag if there have been many failures.
600   continue
      iphase = 1

      ! Restore X, PHI, PSI
      t = told
      if (kp1 >= nsp1) then
         do j = nsp1, kp1
            temp1 = 1/beta(j)
            phi(:, j) = temp1*phi(:, j)
         end do
      end if

630   continue
      psi(1:kp1-1) = psi(2:kp1) - h

      ! Test whether failure is due to nonlinear solver or error test
      if (iernls == 0) goto 660
      iwm(iwloc_ncfn) = iwm(iwloc_ncfn) + 1

      ! The nonlinear solver failed to converge.
      ! Determine the cause of the failure and take appropriate action.
      ! If IERNLS < 0, then return.  Otherwise, reduce the stepsize
      ! and try again, unless too many failures have occurred.
      if (iernls < 0) goto 675
      ncf = ncf + 1
      r = 0.25_rk
      h = h*r
      if (ncf < 10 .and. abs(h) >= hmin) goto 690
      if (idid == 1) idid = -7
      if (nef >= 3) idid = -9
      goto 675

      ! The nonlinear solver converged, and the cause of the failure was the error estimate
      ! exceeding the tolerance.
660   continue
      nef = nef + 1
      iwm(iwloc_netf) = iwm(iwloc_netf) + 1
      if (nef > 1) goto 665

      ! On first error test failure, keep current order or lower order by one.
      ! Compute new stepsize based on differences of the solution.
      k = knew
      temp2 = real(k + 1, rk)
      r = 0.9_rk*(2*est + 1e-4_rk)**(-1/temp2)
      r = max(0.25_rk, min(0.9_rk, r))
      h = h*r
      if (abs(h) >= hmin) goto 690
      idid = -6
      goto 675

      ! On second error test failure, use the current order or decrease order by one.
      ! Reduce the stepsize by a factor of one quarter.
665   continue
      if (nef > 2) goto 670
      k = knew
      r = 0.25_rk
      h = r*h
      if (abs(h) >= hmin) goto 690
      idid = -6
      goto 675

      ! On third and subsequent error test failures, set the order to one,
      ! and reduce the stepsize by a factor of one quarter.
670   continue
      k = 1
      r = 0.25_rk
      h = r*h
      if (abs(h) >= hmin) goto 690
      idid = -6
      goto 675

      ! For all crashes, restore Y to its last value,
      ! interpolate to find YPRIME at last X, and return.

      ! Before returning, verify that the user has not set IDID to a nonnegative value.
      ! If the user has set IDID to a nonnegative value, then reset IDID to be -7,
      ! indicating a failure in the nonlinear system solver.
675   continue
      call ddatrp(t, t, y, ydot, neq, k, phi, psi)
      jstart = 1
      if (idid >= 0) idid = -7
      return

      ! Go back and try this step again.
      ! If this is the first step, reset PSI(1) and rescale PHI(*,2).
690   continue
      if (kold == 0) then
         psi(1) = h
         phi(:, 2) = r*phi(:, 2)
      end if

   end do

end subroutine ddstp

subroutine dcnstr(neq, y, ynew, icnstr, tau, rlx, iret, ivar)
!! This subroutine checks for constraint violations in the proposed new approximate solution
!! `ynew`. If a constraint violation occurs, then a new step length, `tau`, is calculated,
!! and this value is given to the linesearch routine to calculate a new approximate
!! solution `ynew`.

   use daskr_kinds, only: rk, zero
   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(in) :: y(neq)
      !! Current approximate solution vector.
   real(rk), intent(in) :: ynew(neq)
      !! New approximate solution vector.
   integer, intent(in) :: icnstr(neq)
      !! Flags for checking constraints. See [[dcnst0]] for details.
   real(rk), intent(in) :: rlx
      !! Factor for restricting the update size for `icnstr = 2` or `icnstr = -2`.
   real(rk), intent(inout) :: tau
      !! On entry, the current step length for the linesearch. On return, the
      !! adjusted step length if a constraint violation occurred.
   integer, intent(out) :: iret
      !! Output flag.
      !! `0`: `ynew` satisfied all constraints.
      !! `1`: `ynew` did not satisfy all the constraints, requiring a new linesearch.
   integer, intent(out) :: ivar
      !! Index of variable causing the constraint to be violated.

   real(rk), parameter :: fac = 0.6_rk, fac2 = 0.9_rk

   real(rk) :: rdy, rdymx
   integer :: i

   ! Check constraints for proposed new step YNEW.  If a constraint has been violated,
   ! then calculate a new step length, TAU, to be used in the linesearch routine.
   iret = 0
   rdymx = zero
   ivar = 0

   do i = 1, neq

      if (icnstr(i) .eq. 2) then
         rdy = abs((ynew(i) - y(i))/y(i))
         if (rdy .gt. rdymx) then
            rdymx = rdy
            ivar = i
         end if
         if (ynew(i) .le. zero) then
            tau = fac*tau
            ivar = i
            iret = 1
            return
         end if

      elseif (icnstr(i) .eq. 1) then
         if (ynew(i) .lt. zero) then
            tau = fac*tau
            ivar = i
            iret = 1
            return
         end if

      elseif (icnstr(i) .eq. -1) then
         if (ynew(i) .gt. zero) then
            tau = fac*tau
            ivar = i
            iret = 1
            return
         end if

      elseif (icnstr(i) .eq. -2) then
         rdy = abs((ynew(i) - y(i))/y(i))
         if (rdy .gt. rdymx) then
            rdymx = rdy
            ivar = i
         end if
         if (ynew(i) .ge. zero) then
            tau = fac*tau
            ivar = i
            iret = 1
            return
         end if

      end if
   end do

   if (rdymx .ge. rlx) then
      tau = fac2*tau*rlx/rdymx
      iret = 1
   end if

end subroutine dcnstr

subroutine dcnst0(neq, y, icnstr, iret)
!! This routine checks for constraint violations in the initial approximate solution `y`.

   use daskr_kinds, only: rk, zero
   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(in) :: y(neq)
      !! Initial approximate solution.
   integer, intent(in) :: icnstr(neq)
      !! Array indicating which variables are constrained and how. For any element of
      !! `icnstr`, the meaning is as follows. 
      !!  `2`: `y(i) > 0`.
      !!  `1`: `y(i) >= 0`.
      !! `-1`: `y(i) <= 0`.
      !! `-2`: `y(i) < 0`.
      !!  `0`: `y(i)` is unconstrained.
   integer, intent(out) :: iret
      !! Output flag.
      !!     `0`: All constraints satisfied.
      !! `i > 0`: Constraint `i` violated.

   integer :: i

   iret = 0
   do i = 1, neq
      if (icnstr(i) .eq. 2) then
         if (y(i) .le. zero) then
            iret = i
            exit
         end if
      elseif (icnstr(i) .eq. 1) then
         if (y(i) .lt. zero) then
            iret = i
            exit
         end if
      elseif (icnstr(i) .eq. -1) then
         if (y(i) .gt. zero) then
            iret = i
            exit
         end if
      elseif (icnstr(i) .eq. -2) then
         if (y(i) .ge. zero) then
            iret = i
            exit
         end if
      end if
   end do

end subroutine dcnst0

pure subroutine ddawts(neq, iwt, rtol, atol, y, wt, rpar, ipar)
!! This subroutine sets the error weight vector, `wt`, according to:
!!```
!!     wt(i) = rtol(i) * abs(y(i)) + atol(i)
!!```
!! where `rtol` and `atol` are scalars if `iwt = 0`, and vectors if `iwt = 1`.

   use daskr_kinds, only: rk

   integer, intent(in) :: neq
      !! Problem size.
   integer, intent(in) :: iwt
      !! Flag indicating whether `rtol` and `atol` are scalars or vectors.
      !! If `iwt = 0`, then `rtol` and `atol` are scalars.
      !! If `iwt = 1`, then `rtol` and `atol` are vectors.
   real(rk), intent(in) :: rtol(*)
      !! Relative error tolerance.
   real(rk), intent(in) :: atol(*)
      !! Absolute error tolerance.
   real(rk), intent(in) :: y(neq)
      !! Solution vector.
   real(rk), intent(out) :: wt(neq)
      !! Error weight vector.
   real(rk), intent(in) :: rpar(*)
      !! User real workspace.
   integer, intent(in) :: ipar(*)
      !! User integer workspace.

   if (iwt == 0) then
      wt = rtol(1)*abs(y) + atol(1)
   else
      wt = rtol(1:neq)*abs(y) + atol(1:neq)
   end if

end subroutine ddawts

pure subroutine dinvwt(neq, wt, ierr)
!! This subroutine checks the error weight vector `wt` for components that are less
!! or equal to 0, and if none are found, it inverts `wt` in place. This replaces
!! division operations with multiplications in all norm evaluations.

   use daskr_kinds, only: rk, zero
   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(inout) :: wt(neq)
      !! Weights for error control.
   integer, intent(out) :: ierr
      !! Completion flag.
      !! `0`: all elements were found positive.
      !! `i > 0`: `w(i)` is the first non positive weight.

   integer :: i

   do i = 1, neq
      if (wt(i) .le. zero) then
         ierr = i
         return
      end if
   end do

   wt = 1/wt
   ierr = 0

end subroutine dinvwt

pure subroutine ddatrp(t, tout, yout, ydotout, neq, kold, phi, psi)
!! The methods in subroutine [[ddstp]] use polynomials to approximate the solution.
!! This routine approximates the solution and its derivative at `tout` by evaluating
!! one of these polynomials, and its derivative, there. Information defining the
!! polynomial is passed from the caller.

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

   real(rk), intent(in) :: t
      !! Current value of the independent variable.
   real(rk), intent(in) :: tout
      !! Value of the independent variable at which the solution is desired.
   real(rk), intent(out) :: yout(neq)
      !! Interpolated approximation to `y` at `tout`.
   real(rk), intent(out) :: ydotout(neq)
      !! Interpolated approximation to `ydot` at `tout`.
   integer, intent(in) :: neq
      !! Problem size.
   integer, intent(in) :: kold
      !! Order used on last successful step.
   real(rk), intent(in) :: phi(neq, kold + 1)
      !! Array of scaled divided differences of `y`.
   real(rk), intent(in) :: psi(kold + 1)
      !! Array of past stepsize history.

   integer :: j
   real(rk) :: c, d, gama, dt

   dt = tout - t
   yout = phi(:, 1)
   ydotout = zero
   c = one
   d = zero
   gama = dt/psi(1)

   do j = 2, kold + 1
      d = d*gama + c/psi(j - 1)
      c = c*gama
      gama = (dt + psi(j - 1))/psi(j)
      yout = yout + c*phi(:, j)
      ydotout = ydotout + d*phi(:, j)
   end do

end subroutine ddatrp

real(rk) pure function ddwnrm(neq, v, rwt, rpar, ipar)
!! This function computes the weighted root-mean-square norm of the vector `v` with
!! reciprocal weights `rwt`:
!!
!! $$ \sqrt{ \frac{1}{N} \sum_{i=1}^N \left( v_i \cdot w_i \right)^2 } $$

   use daskr_kinds, only: rk, zero
   implicit none

   integer, intent(in) :: neq
      !! Vector length.
   real(rk), intent(in) :: v(neq)
      !! Vector.
   real(rk), intent(in) :: rwt(neq)
      !! Reciprocal weights.
   real(rk), intent(in) :: rpar(*)
      !! User real workspace.
   integer, intent(in) :: ipar(*)
      !! User integer workspace.

   integer :: i
   real(rk) :: accum, vmax

   ddwnrm = zero
   vmax = zero
   do i = 1, neq
      if (abs(v(i)*rwt(i)) .gt. vmax) then
         vmax = abs(v(i)*rwt(i))
      end if
   end do

   if (vmax .le. zero) return

   accum = zero
   do i = 1, neq
      accum = accum + ((v(i)*rwt(i))/vmax)**2
   end do

   ddwnrm = vmax*sqrt(accum/neq)

end function ddwnrm

subroutine ddasid( &
   t, y, ydot, neq, icopt, idy, res, jac, dum1, h, tscale, &
   wt, dum2, rpar, ipar, dum3, delta, r, y0, ydot0, dum4, rwm, iwm, cj, uround, &
   dum5, dum6, dum7, epscon, ratemax, stptol, dum8, icnflg, icnstr, iernls)
!! This routine solves a nonlinear system of algebraic equations of the form:
!!
!!  $$ G(t, y, \dot{y}) = 0 $$
!!
!! for the unknown parts of \(y\) and \(\dot{y}\) in the initial conditions. 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
   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.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   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
      !! Scaling factor for this initial condition calculation.
   real(rk), intent(in) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: wt(neq)
      !! Weights for error criterion.
   integer, intent(inout) :: dum2
      !! Dummy argument.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(out) :: dum3(neq)
      !! Dummy argument.
   real(rk), intent(inout) :: delta(neq)
      !! Real work array.
   real(rk), intent(out) :: r(neq)
      !! Real work array.
   real(rk), intent(out) :: y0(neq)
      !! Real work array.
   real(rk), intent(out) :: ydot0(neq)
      !! Real work array.
   real(rk), intent(inout) :: dum4(neq)
      !! Dummy argument.
   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) :: uround ! @todo: used?
      !! Unit roundoff.
   real(rk), intent(in) :: dum5
      !! Dummy argument.
   real(rk), intent(in) :: dum6
      !! Dummy argument.
   real(rk), intent(in) :: dum7
      !! Dummy argument.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   real(rk), intent(in) :: ratemax
      !! Maximum convergence rate for which Newton iteration is considered converging.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(in) :: dum8
      !! Dummy argument.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   integer, intent(out) :: iernls
      !! Error flag for nonlinear solver.
      !! `0`: nonlinear solver converged.
      !! `1`, `2`: recoverable error inside nonlinear solver.
      !! `1`: retry with current (y,y').
      !! `2`: retry with original (y,y').
      !! `-1`: unrecoverable error in nonlinear solver.

   integer :: ierj, iernew, ires, mxnit, mxnj, nj
   logical :: failed

   ! Perform initializations.
   mxnit = iwm(iwloc_mxnit)
   mxnj = iwm(iwloc_mxnj)
   iernls = 0
   nj = 0
   failed = .false.

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

   if (ires .lt. 0) then
      failed = .true.
   end if

   ! Looping point for updating the Jacobian.
   if (.not. failed) then
      do

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

         ! Reevaluate the iteration matrix, J = dG/dY + CJ*dG/dYDOT.
         call dmatd(neq, t, y, ydot, delta, cj, h, ierj, wt, r, rwm, iwm, res, ires, &
                    uround, jac, rpar, ipar)
         nj = nj + 1
         iwm(iwloc_nje) = iwm(iwloc_nje) + 1

         if ((ires .lt. 0) .or. (ierj .ne. 0)) then
            failed = .true.
            exit
         end if

         ! Call the nonlinear Newton solver.
         call dnsid(t, y, ydot, neq, icopt, idy, res, wt, rpar, ipar, delta, r, &
                    y0, ydot0, rwm, iwm, cj, tscale, epscon, ratemax, mxnit, stptol, &
                    icnflg, icnstr, iernew)

         if ((iernew .eq. 1) .and. (nj .lt. mxnj)) then
            ! MXNIT iterations were done, the convergence rate is < 1,
            ! and the number of Jacobian evaluations is less than MXNJ.
            ! Call RES, reevaluate the Jacobian, and try again.
            ires = 0
            call res(t, y, ydot, cj, delta, ires, rpar, ipar)
            iwm(iwloc_nre) = iwm(iwloc_nre) + 1
            if (ires .lt. 0) then
               failed = .true.
               exit
            end if
            cycle
         else
            exit ! Either success or convergence failure
         end if

      end do
   end if

   ! Unsuccessful exits from nonlinear solver.
   ! Compute IERNLS accordingly.
   if (failed) then
      iernls = 2
      if (ires .le. -2) iernls = -1
      return
   end if

   if (iernew .ne. 0) then
      iernls = min(iernew, 2)
   end if

end subroutine ddasid

subroutine dnsid( &
   t, y, ydot, neq, icopt, idy, res, wt, rpar, ipar, &
   delta, r, y0, ydot0, rwm, iwm, cj, tscale, &
   epscon, ratemax, maxit, stptol, icnflg, icnstr, iernew)
!! This routine solves a nonlinear system of algebraic equations of the form:
!!
!! $$ G(t, y, \dot{y}) = 0 $$
!!
!! for the unknown parts of \(y\) and \(\dot{y}\) in the initial conditions. 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.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   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(inout) :: delta(neq)
      !! Residual vector on entry, and real workspace.
   real(rk), intent(out) :: r(neq)
      !! Real work array.
   real(rk), intent(out) :: y0(neq)
      !! Real work array.
   real(rk), intent(out) :: ydot0(neq)
      !! Real work array.
   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) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   real(rk), intent(in) :: ratemax
      !! Maximum convergence rate for which Newton iteration is considered converging.
   integer, intent(in) :: maxit
      !! Maximum number of Newton iterations allowed.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   integer, intent(out) :: iernew
      !! Error flag for the Newton iteration.
      !!  `0`: the Newton iteration converged.
      !!  `1`: failed to converge, but `rate < ratemax`.
      !!  `2`: failed to converge, but `rate > ratemax`.
      !!  `3`: other recoverable error (`ires = -1` or linesearch failed).
      !! `-1`: unrecoverable error inside Newton iteration.

   integer :: ires, iret, lsoff, m
   real(rk) :: deltanorm, fnorm, oldfnorm, rate, rlx
   real(rk), external :: ddwnrm

   ! Initializations.
   lsoff = iwm(iwloc_lsoff)
   rate = one
   rlx = 0.4_rk
   iernew = 0

   ! Compute a new step vector DELTA by back-substitution.
   call dslvd(neq, delta, rwm, iwm)

   ! Get norm of DELTA. Return now if norm(DELTA) .le. EPCON.
   deltanorm = ddwnrm(neq, delta, wt, rpar, ipar)
   fnorm = deltanorm
   if (tscale .gt. zero) fnorm = fnorm*tscale*abs(cj)
   if (fnorm .le. epscon) return

   ! Newton iteration loop.
   m = 0
   do

      iwm(iwloc_nni) = iwm(iwloc_nni) + 1

      ! Call linesearch routine for global strategy and set RATE
      oldfnorm = fnorm
      call dlinsd(neq, y, t, ydot, cj, tscale, delta, deltanorm, wt, &
                  lsoff, stptol, iret, res, ires, rwm, iwm, fnorm, icopt, &
                  idy, r, y0, ydot0, icnflg, icnstr, rlx, rpar, ipar)
      rate = fnorm/oldfnorm

      ! Check for error condition from linesearch.
      if (iret .ne. 0) goto 390

      ! Test for convergence of the iteration, and return or loop.
      if (fnorm .le. epscon) return

      ! The iteration has not yet converged. Update M.
      ! Test whether the maximum number of iterations have been tried.
      m = m + 1
      if (m .ge. maxit) goto 380

      ! Copy the residual to DELTA and its norm to DELNRM, and loop for another iteration.
      delta = r
      deltanorm = fnorm

   end do

   ! The maximum number of iterations was done. Set IERNEW and return.
380 continue
   if (rate .le. ratemax) then
      iernew = 1
   else
      iernew = 2
   end if
   return

   ! Errors
390 continue
   if (ires .le. -2) then
      iernew = -1
   else
      iernew = 3
   end if

end subroutine dnsid

subroutine dlinsd( &
   neq, y, t, ydot, cj, tscale, p, pnorm, wt, &
   lsoff, stptol, iret, res, ires, rwm, iwm, &
   fnorm, icopt, idy, r, ynew, ydotnew, icnflg, &
   icnstr, rlx, rpar, ipar)
!! This routine uses a linesearch algorithm to calculate a new \( (y,\dot{y}) \) pair,
!! denoted \( (y_{new},\dot{y}_{new}) \), such that:
!!
!!   $$ f(y_{new},\dot{y}_{new}) \leq (1 - 2 \alpha \lambda) f(y,\dot{y})$$
!!
!! where \(0 < \lambda \le 1\), and \( f(y,\dot{y}) \) is defined as:
!!
!! $$ f(y,\dot{y}) = (1/2) \| (J^{-1} G(t,y,\dot{y}) \|^2 $$
!!
!! where the norm is the weighted root mean square vector norm, \(G\) is the DAE system
!! residual function, and \(J\) is the system iteration matrix (Jacobian).

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

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(inout) :: y(neq)
      !! On entry, the current solution vector. On return, the new solution vector.
   real(rk), intent(in) :: t
      !! Independent variable.
   real(rk), intent(inout) :: ydot(neq)
      !! On entry, the current derivative vector. On return, the new derivative vector.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: p(neq)
      !! Approximate Newton step used in backtracking.
   real(rk), intent(inout) :: pnorm
      !! Weighted root mean square norm of `p`.
   real(rk), intent(inout) :: wt(neq)
      !! Scaling factors.
   integer, intent(in) :: lsoff ! @todo: convert to logical
      !! Flag showing whether the linesearch algorithm is to be invoked.
      !! `0`: do the linesearch.
      !! `1`: turn off linesearch.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(out) :: iret
      !! Return flag.
      !! `0`: a satisfactory (y,y') was found.
      !! `1`: the routine failed to find a new (y,y') pair that was sufficiently distinct
      !! from the current one.
      !! `2`: a failure in `res`.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   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(inout) :: fnorm
      !! Value of \( \sqrt{2f} \) for the current \( (y,\dot{y}) \) on input and output.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   real(rk), intent(out) :: r(neq)
      !! Scaled residual vector \(J^{-1} G(t,y,\dot{y})\).
   real(rk), intent(out) :: ynew(neq)
      !! New solution vector.
   real(rk), intent(out) :: ydotnew(neq)
      !! New derivative of solution vector.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   real(rk), intent(in) :: rlx
      !! Factor for restricting update size in [[dcnstr]].
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   real(rk), parameter :: alpha = 1e-4_rk
   integer :: ivar, kprint
   real(rk) :: f1norm, f1normp, fnormp, ratio, ratio1, rl, rlmin, slpi, tau
   character(len=80) :: msg

   kprint = iwm(iwloc_kprint)
   f1norm = (fnorm**2)/2
   ratio = one

   if (kprint .ge. 2) then
      msg = '------ IN ROUTINE DLINSD-- PNRM = (R1)'
      call xerrwd(msg, 38, 901, 0, 0, 0, 0, 1, pnorm, zero)
   end if
   tau = pnorm
   rl = one

   ! Check for violations of the constraints, if any are imposed.
   ! If any violations are found, the step vector P is rescaled, and the
   ! constraint check is repeated, until no violations are found.
   if (icnflg .ne. 0) then
      do
         call dyypnw(neq, y, ydot, cj, rl, p, icopt, idy, ynew, ydotnew)
         call dcnstr(neq, y, ynew, icnstr, tau, rlx, iret, ivar)
         if (iret /= 1) exit
         ratio1 = tau/pnorm
         ratio = ratio*ratio1
         p = p*ratio1
         pnorm = tau
         if (kprint .ge. 2) then
            msg = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)'
            call xerrwd(msg, 50, 902, 0, 1, ivar, 0, 1, pnorm, zero)
         end if
         if (pnorm .le. stptol) then
            iret = 1
            return
         end if
      end do
   end if

   slpi = -2*f1norm*ratio
   rlmin = stptol/pnorm
   if ((lsoff .eq. 0) .and. (kprint .ge. 2)) then
      msg = '------ MIN. LAMBDA = (R1)'
      call xerrwd(msg, 25, 903, 0, 0, 0, 0, 1, rlmin, zero)
   end if

   ! Begin iteration to find RL value satisfying alpha-condition.
   ! If RL becomes less than RLMIN, then terminate with IRET = 1.
   do
      call dyypnw(neq, y, ydot, cj, rl, p, icopt, idy, ynew, ydotnew)
      call dfnrmd(neq, ynew, t, ydotnew, r, cj, tscale, wt, res, ires, &
                  fnormp, rwm, iwm, rpar, ipar)
      iwm(iwloc_nre) = iwm(iwloc_nre) + 1
      if (ires .ne. 0) then
         iret = 2
         return
      end if

      if (lsoff .eq. 1) goto 150

      f1normp = (fnormp**2)/2
      if (kprint .ge. 2) then
         msg = '------ LAMBDA = (R1)'
         call xerrwd(msg, 20, 904, 0, 0, 0, 0, 1, rl, zero)
         msg = '------ NORM(F1) = (R1),  NORM(F1NEW) = (R2)'
         call xerrwd(msg, 43, 905, 0, 0, 0, 0, 2, f1norm, f1normp)
      end if

      if (f1normp .gt. f1norm + alpha*slpi*rl) goto 200

      ! Alpha-condition is satisfied, or linesearch is turned off.
      ! Copy YNEW,YDOTNEW to Y,YDOT and return.
150   continue
      iret = 0
      y = ynew
      ydot = ydotnew
      fnorm = fnormp
      if (kprint .ge. 1) then
         msg = '------ LEAVING ROUTINE DLINSD, FNRM = (R1)'
         call xerrwd(msg, 42, 906, 0, 0, 0, 0, 1, fnorm, zero)
      end if
      return

      ! Alpha-condition not satisfied. Perform backtrack to compute new RL
      ! value. If no satisfactory YNEW,YDOTNEW can be found sufficiently
      ! distinct from Y,YDOT, then return IRET = 1.
200   continue
      if (rl .lt. rlmin) then
         iret = 1
         return
      end if

      rl = rl/2

   end do

end subroutine dlinsd

subroutine dfnrmd( &
   neq, y, t, ydot, r, cj, tscale, wt, &
   res, ires, rnorm, rwm, iwm, rpar, ipar)
!! This routine calculates the scaled preconditioned norm of the nonlinear function used
!! in the nonlinear iteration for obtaining consistent initial conditions. Specifically,
!! it calculates the weighted root-mean-square norm of the vector:
!!
!! $$  r = J^{-1} G(t,y,\dot{y}) $$
!!
!! where \(J\) is the Jacobian matrix and \(G\) is the DAE equation vector.

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

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(in) :: y(neq)
      !! Solution vector.
   real(rk), intent(in) :: t
      !! Independent variable.
   real(rk), intent(in) :: ydot(neq)
      !! Derivative of solution vector.
   real(rk), intent(out) :: r(neq)
      !! Result vector.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: tscale ! @todo: what is "t?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: wt(neq)
      !! Scaling factors.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   real(rk), intent(out) :: rnorm
      !! Weighted root-mean-square norm of `r`.
   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(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   real(rk), external :: ddwnrm ! @todo: remove this once inside module

   ! Call RES routine.
   ires = 0
   call res(t, y, ydot, cj, r, ires, rpar, ipar)
   if (ires .lt. 0) return

   ! Apply inverse of Jacobian to vector R.
   call dslvd(neq, r, rwm, iwm)

   ! Calculate norm of R.
   rnorm = ddwnrm(neq, r, wt, rpar, ipar)
   if (tscale .gt. zero) rnorm = rnorm*tscale*abs(cj)

end subroutine dfnrmd

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

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

subroutine dmatd( &
   neq, t, y, ydot, delta, cj, h, ierr, ewt, e, &
   rwm, iwm, res, ires, uround, jac, rpar, ipar)
!! This routine computes the iteration matrix:
!!
!! $$ J = \partial G/ \partial y + c_J \partial G/ \partial \dot{y} $$
!!
!! Here \(J\) is computed by the user-supplied routine `jac` if `mtype` is 1 or 4, or
!! by numerical difference quotients if `mtype` is 2 or 5.

   use daskr_kinds, only: rk, zero
   use daskr
   use linpack, only: dgefa, dgbfa
   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   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.
   real(rk), intent(inout) :: delta(neq)
      !! Residual evaluated at (t,y,y').
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: h
      !! Current stepsize in integration.
   integer, intent(out) :: ierr
      !! Error flag returned by [[dgefa]] and [[dgbfa]].
      !!     `0`: matrix is not singular.
      !! `k > 0`: matrix is singular.
   real(rk), intent(inout) :: ewt(neq)
      !! Vector of error weights for computing norms.
   real(rk), intent(inout) :: e(neq)
      !! Real work array.
   real(rk), intent(inout) :: rwm(*)
      !! Real workspace for matrices. On output, it contains the LU decomposition
      !! of the iteration matrix.
   integer, intent(inout) :: iwm(*)
      !! Integer workspace containing matrix information.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   real(rk), intent(in) :: uround ! @todo: remove this
      !! Unit roundoff error.
   procedure(jac_t) :: jac
      !! User-defined Jacobian routine.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   integer :: i, i1, i2, ii, ipsave, isave, j, k, l, lenpd, lipvt, mba, mband, meb1, &
              meband, msave, mtype, n, nrow
   real(rk) :: del, delinv, squr, ypsave, ysave
   logical :: isdense

   lipvt = iwm(iwloc_liwp)
   mtype = iwm(iwloc_mtype)
   ierr = 0

   select case (mtype)
   case (1)
      ! Dense user-supplied matrix.
      isdense = .true.
      lenpd = iwm(iwloc_npd)
      rwm(1:lenpd) = zero
      call jac(t, y, ydot, rwm, cj, rpar, ipar)

   case (2)
      ! Dense finite-difference-generated matrix.
      isdense = .true.
      ires = 0
      nrow = 0
      squr = sqrt(uround)
      do i = 1, neq
         del = max(squr*max(abs(y(i)), abs(h*ydot(i))), 1/ewt(i))
         del = sign(del, h*ydot(i))
         del = (y(i) + del) - y(i)
         ysave = y(i)
         ypsave = ydot(i)
         y(i) = y(i) + del
         ydot(i) = ydot(i) + cj*del
         call res(t, y, ydot, cj, e, ires, rpar, ipar)
         iwm(iwloc_nre) = iwm(iwloc_nre) + 1
         if (ires .lt. 0) return
         delinv = 1/del
         do l = 1, neq
            rwm(nrow + l) = (e(l) - delta(l))*delinv
         end do
         nrow = nrow + neq
         y(i) = ysave
         ydot(i) = ypsave
      end do

   case (3)
      ! dummy section for mtype=3.
      error stop "error: mtype=3 not implemented"

   case (4)
      ! Banded user-supplied matrix.
      isdense = .false.
      lenpd = iwm(iwloc_npd)
      rwm(1:lenpd) = zero
      call jac(t, y, ydot, rwm, cj, rpar, ipar)
      meband = 2*iwm(iwloc_ml) + iwm(iwloc_mu) + 1

   case (5)
      ! Banded finite-difference-generated matrix.
      isdense = .false.
      mband = iwm(iwloc_ml) + iwm(iwloc_mu) + 1
      mba = min(mband, neq)
      meband = mband + iwm(iwloc_ml)
      meb1 = meband - 1
      msave = (neq/mband) + 1
      isave = iwm(iwloc_npd)
      ipsave = isave + msave
      ires = 0
      squr = sqrt(uround)
      do j = 1, mba
         do n = j, neq, mband
            k = (n - j)/mband + 1
            rwm(isave + k) = y(n)
            rwm(ipsave + k) = ydot(n)
            del = max(squr*max(abs(y(n)), abs(h*ydot(n))), 1/ewt(n))
            del = sign(del, h*ydot(n))
            del = (y(n) + del) - y(n)
            y(n) = y(n) + del
            ydot(n) = ydot(n) + cj*del
         end do
         call res(t, y, ydot, cj, e, ires, rpar, ipar)
         iwm(iwloc_nre) = iwm(iwloc_nre) + 1
         if (ires .lt. 0) return
         do n = j, neq, mband
            k = (n - j)/mband + 1
            y(n) = rwm(isave + k)
            ydot(n) = rwm(ipsave + k)
            del = max(squr*max(abs(y(n)), abs(h*ydot(n))), 1/ewt(n))
            del = sign(del, h*ydot(n))
            del = (y(n) + del) - y(n)
            delinv = 1/del
            i1 = max(1, (n - iwm(iwloc_mu)))
            i2 = min(neq, (n + iwm(iwloc_ml)))
            ii = n*meb1 - iwm(iwloc_ml)
            do i = i1, i2
               rwm(ii + i) = (e(i) - delta(i))*delinv
            end do
         end do
      end do

   case default
      error stop "error: unexpected value for mtype"

   end select

   if (isdense) then
      call dgefa(rwm, neq, neq, iwm(lipvt), ierr)
   else
      call dgbfa(rwm, meband, neq, iwm(iwloc_ml), iwm(iwloc_mu), iwm(lipvt), ierr)
   end if

end

subroutine dslvd(neq, delta, rwm, iwm)
!! This routine manages the solution of the linear system arising in the Newton iteration.
!! Real matrix information and real temporary storage is stored in the array `rwm`, while
!! integer matrix information is stored in the array `iwm`.
!! For a dense matrix, the LINPACK routine [[dgesl]] is called. For a banded matrix, the
!! LINPACK routine [[dgbsl]] is called.

   use daskr_kinds, only: rk
   use daskr
   use linpack, only: dgesl, dgbsl
   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(inout) :: delta(neq)
      !! On entry, the right-hand side vector of the linear system to be solved,
      !! and, on exit, the solution vector.
   real(rk), intent(inout) :: rwm(*)
      !! Real workspace for the linear system solver.
   integer, intent(inout) :: iwm(*)
      !! Integer workspace for the linear system solver.

   integer :: lipvt, meband, mtype

   lipvt = iwm(iwloc_liwp)
   mtype = iwm(iwloc_mtype)

   select case (mtype)
   case (1, 2)
      ! dense matrix.
      call dgesl(rwm, neq, neq, iwm(lipvt), delta, 0)
   case (3)
      ! dummy section for mtype=3.
      error stop "error: mtype=3 not implemented"
   case (4, 5)
      ! banded matrix.
      meband = 2*iwm(iwloc_ml) + iwm(iwloc_mu) + 1
      call dgbsl(rwm, meband, neq, iwm(iwloc_ml), iwm(iwloc_mu), iwm(lipvt), delta, 0)
   case default
      error stop "error: unexpected value for mtype"
   end select

end subroutine dslvd

subroutine ddasik( &
   t, y, ydot, neq, icopt, idy, res, jack, psol, h, tscale, &
   wt, jskip, rpar, ipar, savr, delta, r, y0, ydot0, pwk, rwm, iwm, cj, uround, &
   epli, sqrtn, rsqrtn, epscon, ratemax, stptol, jflg, icnflg, icnstr, iernls)
!! This routine solves a nonlinear system of algebraic equations of the form:
!!
!!  $$ G(t, y, \dot{y}) = 0 $$
!!
!! for the unknown parts of \(y\) and \(\dot{y}\) in the initial conditions. An initial
!! value for \(y\) and initial guess for \(\dot{y}\) are input. The method used is a
!! Newton scheme with Krylov iteration and a linesearch algorithm.

   use daskr_kinds, only: rk
   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.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(jack_t) :: jack
      !! User-supplied routine to update the preconditioner.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   real(rk), intent(in) :: h
      !! Scaling factor for this initial condition calculation.
   real(rk), intent(in) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: wt(neq)
      !! Weights for error criterion.
   integer, intent(inout) :: jskip
      !! Flag to signal if initial `jack` call is to be skipped.
      !! `1`: skip the call,
      !! `0`: do not skip call.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(out) :: savr(neq)
      !! Real work array.
   real(rk), intent(inout) :: delta(neq)
      !! Real work array.
   real(rk), intent(out) :: r(neq)
      !! Real work array.
   real(rk), intent(out) :: y0(neq)
      !! Real work array.
   real(rk), intent(out) :: ydot0(neq)
      !! Real work array.
   real(rk), intent(inout) :: pwk(neq)
      !! Real work array.
   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) :: uround
      !! Unit roundoff (not used).
   real(rk), intent(in) :: epli ! @note: not the same as 'epslin'!
      !! Convergence test constant for the iterative solution of linear systems. See
      !! [[daskr]] for details.
   real(rk), intent(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   real(rk), intent(in) :: ratemax
      !! Maximum convergence rate for which Newton iteration is considered converging.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(in) :: jflg
      !! Flag indicating whether a `jac` routine is supplied by the user.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   integer, intent(out) :: iernls
      !! Error flag for nonlinear solver.
      !! `0`: nonlinear solver converged.
      !! `1`, `2`: recoverable error inside nonlinear solver.
      !! `1`: retry with current (y,y').
      !! `2`: retry with original (y,y').
      !! `-1`: unrecoverable error in nonlinear solver.

   integer :: iernew, ierpj, ires, liwp, lrwp, maxnit, maxnj, nj
   real(rk) :: epslin
   logical :: failed

   ! Perform initializations.
   lrwp = iwm(iwloc_lrwp)
   liwp = iwm(iwloc_liwp)
   maxnit = iwm(iwloc_mxnit)
   maxnj = iwm(iwloc_mxnj)
   nj = 0
   epslin = epli*epscon
   failed = .false.

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

   ! Preconditioner update loop
   if (.not. failed) then
      do

         ! Set all error flags to zero.
         iernls = 0
         ierpj = 0
         ires = 0
         iernew = 0

         ! If a Jacobian routine was supplied, call it.
         if ((jflg .eq. 1) .and. (jskip .eq. 0)) then
            nj = nj + 1
            iwm(iwloc_nje) = iwm(iwloc_nje) + 1
            call jack(res, ires, neq, t, y, ydot, wt, delta, r, h, cj, &
                      rwm(lrwp), iwm(liwp), ierpj, rpar, ipar)

            if ((ires .lt. 0) .or. (ierpj .ne. 0)) then
               failed = .true.
               exit
            end if
         end if

         jskip = 0

         ! Call the nonlinear Newton solver.
         call dnsik(t, y, ydot, neq, icopt, idy, res, psol, wt, rpar, ipar, &
                    savr, delta, r, y0, ydot0, pwk, rwm, iwm, cj, tscale, sqrtn, rsqrtn, &
                    epslin, epscon, ratemax, maxnit, stptol, icnflg, icnstr, iernew)

         if ((iernew .eq. 1) .and. (nj .lt. maxnj) .and. (jflg .eq. 1)) then
            ! Up to MXNIT iterations were done, the convergence rate is < 1,
            ! a Jacobian routine is supplied, and the number of JACK calls is less than MXNJ.
            ! Copy the residual SAVR to DELTA, call JACK, and try again.
            delta = savr
            cycle
         else
            exit ! Either success or convergence failure
         end if

      end do
   end if

   if (failed) then
      iernls = 2
      if (ires .le. -2) iernls = -1
      return
   end if

   if (iernew .ne. 0) then
      iernls = min(iernew, 2)
   end if

end subroutine ddasik

subroutine dnsik( &
   t, y, ydot, neq, icopt, idy, res, psol, wt, rpar, ipar, &
   savr, delta, r, yic, ydotic, pwk, rwm, iwm, cj, tscale, sqrtn, rsqrtn, epslin, &
   epscon, ratemax, maxit, stptol, icnflg, icnstr, iernew)
!! This routine solves a nonlinear system of algebraic equations of the form:
!!
!!  $$ G(t, y, \dot{y}) = 0 $$
!!
!! for the unknown parts of \(y\) and \(\dot{y}\) in the initial conditions. The method
!! used is a Newton scheme combined with a linesearch algorithm, using Krylov iterative
!! linear system methods.

   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.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   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(out) :: savr(neq)
      !! Real work array.
   real(rk), intent(inout) :: delta(neq)
      !! Residual vector on entry, and real workspace.
   real(rk), intent(out) :: r(neq)
      !! Real work array.
   real(rk), intent(out) :: yic(neq)
      !! Real work array.
   real(rk), intent(out) :: ydotic(neq)
      !! Real work array.
   real(rk), intent(inout) :: pwk(neq)
      !! Real work array.
   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) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   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.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   real(rk), intent(in) :: ratemax
      !! Maximum convergence rate for which Newton iteration is considered converging.
   integer, intent(in) :: maxit
      !! Maximum number of Newton iterations allowed.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   integer, intent(out) :: iernew
      !! Error flag for the Newton iteration.
      !!  `0`: the Newton iteration converged.
      !!  `1`: failed to converge, but `rate < 1`, or the residual norm was reduced by a factor of 0.1.
      !!  `2`: failed to converge, but `rate > ratemx`.
      !!  `3`: other recoverable error.
      !! `-1`: unrecoverable error inside Newton iteration.

   integer :: ierr, iersl, ires, iret, liwp, lsoff, lrwp, m
   real(rk) :: deltanorm, fnorm, fnorm0, oldfnorm, rate, rhok, rlx
   logical :: iserror, ismaxit
   real(rk), external :: ddwnrm

   ! Initializations.
   lsoff = iwm(iwloc_lsoff)
   lrwp = iwm(iwloc_lrwp)
   liwp = iwm(iwloc_liwp)
   rate = one
   rlx = 0.4_rk
   iernew = 0

   ! Save residual in SAVR.
   savr = delta

   ! Compute norm of (P-inverse)*(residual).
   call dfnrmk(neq, y, t, ydot, savr, r, cj, tscale, wt, &
               sqrtn, rsqrtn, res, ires, psol, 1, ierr, fnorm, epslin, &
               rwm(lrwp), iwm(liwp), pwk, rpar, ipar)
   iwm(iwloc_nps) = iwm(iwloc_nps) + 1
   if (ierr .ne. 0) then
      iernew = 3
      return
   end if

   ! Return now if residual norm is .le. EPSCON.
   if (fnorm .le. epscon) return

   ! Newton iteration loop. M is the Newton iteration counter.
   fnorm0 = fnorm
   iserror = .false.
   ismaxit = .false.
   m = 0
   do

      iwm(iwloc_nni) = iwm(iwloc_nni) + 1

      ! Compute a new step vector 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 .ne. 0) .or. (iersl .ne. 0)) then
         iserror = .true.
         exit
      end if

      ! Get norm of DELTA. Return now if DELTA is zero.
      deltanorm = ddwnrm(neq, delta, wt, rpar, ipar)
      if (deltanorm .eq. zero) exit

      ! Call linesearch routine for global strategy and set RATE.
      oldfnorm = fnorm
      call dlinsk(neq, y, t, ydot, savr, cj, tscale, delta, deltanorm, &
                  wt, sqrtn, rsqrtn, lsoff, stptol, iret, res, ires, psol, &
                  rwm, iwm, rhok, fnorm, icopt, idy, rwm(lrwp), iwm(liwp), r, epslin, &
                  yic, ydotic, pwk, icnflg, icnstr, rlx, rpar, ipar)
      rate = fnorm/oldfnorm

      ! Check for error condition from linesearch.
      if (iret .ne. 0) then
         iserror = .true.
         exit
      end if

      ! Test for convergence of the iteration, and return or loop.
      if (fnorm .le. epscon) exit

      ! The iteration has not yet converged. Update M.
      ! Test whether the maximum number of iterations have been tried.
      m = m + 1
      if (m .ge. maxit) then
         ismaxit = .true.
         exit
      end if

      ! Copy the residual SAVR to DELTA and loop for another iteration.
      delta = savr

   end do

   ! The maximum number of iterations was done. Set IERNEW and return.
   if (ismaxit) then
      if ((rate .le. ratemax) .or. (fnorm .le. fnorm0/10)) then
         iernew = 1
      else
         iernew = 2
      end if
   end if

   ! Errors
   if (iserror) then
      if ((ires .le. -2) .or. (iersl .lt. 0)) then
         iernew = -1
      else
         iernew = 3
         if ((ires .eq. 0) .and. (iersl .eq. 1) .and. (m .ge. 2) .and. (rate .lt. one)) then
            iernew = 1
         end if
      end if
   end if

end subroutine dnsik

subroutine dlinsk( &
   neq, y, t, ydot, savr, cj, tscale, p, pnorm, wght, &
   sqrtn, rsqrtn, lsoff, stptol, iret, res, ires, psol, &
   rwm, iwm, rhok, fnorm, icopt, idy, rwp, iwp, r, epslin, ynew, ydotnew, &
   pwk, icnflg, icnstr, rlx, rpar, ipar)
!! This routine uses a linesearch algorithm to calculate a new \( (y,\dot{y}) \) pair,
!! denoted \( (y_{new},\dot{y}_{new}) \), such that:
!!
!!   $$ f(y_{new},\dot{y}_{new}) \leq (1 - 2 \alpha \lambda) f(y,\dot{y})$$
!!
!! where \(0 < \lambda \le 1\), and \( f(y,\dot{y}) \) is defined as:
!!
!! $$ f(y,\dot{y}) = (1/2) \| (P^{-1} G(t,y,\dot{y}) \|^2 $$
!!
!! where the norm is the weighted root mean square vector norm, \(G\) is the DAE system
!! residual function, and \(P\) is the preconditioner used in the Krylov iteration.

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

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(inout) :: y(neq)
      !! On entry, the current solution vector. On return, the new solution vector.
   real(rk), intent(in) :: t
      !! Independent variable.
   real(rk), intent(inout) :: ydot(neq)
      !! On entry, the current derivative vector. On return, the new derivative vector.
   real(rk), intent(out) :: savr(neq)
      !! New residual vector.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: p(neq)
      !! Approximate Newton step used in backtracking.
   real(rk), intent(inout) :: pnorm
      !! Weighted root mean square norm of `p`.
   real(rk), intent(inout) :: wght(neq)
      !! Scaling factors.
   real(rk), intent(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   integer, intent(in) :: lsoff ! @todo: convert to logical
      !! Flag showing whether the linesearch algorithm is to be invoked.
      !! `0`: do the linesearch.
      !! `1`: turn off linesearch.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(out) :: iret
      !! Return flag.
      !! `0`: a satisfactory (y,y') was found.
      !! `1`: the routine failed to find a new (y,y') pair that was sufficiently distinct
      !! from the current one.
      !! `2`: a failure in `res` or `psol`.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   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) :: rhok
      !! Weighted norm of preconditioned Krylov residual (not used).
   real(rk), intent(inout) :: fnorm
      !! Value of \( \sqrt{2f} \) for the current \( (y,\dot{y}) \) on input and output.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   real(rk), intent(inout) :: rwp(*)
      !! Real workspace for the linear system solver.
   integer, intent(inout) :: iwp(*)
      !! Integer workspace for the linear system solver.
   real(rk), intent(out) :: r(neq)
      !! Scaled residual vector \(P^{-1} G(t,y,\dot{y})\).
   real(rk), intent(in) :: epslin
      !! Tolerance for linear system.
   real(rk), intent(out) :: ynew(neq)
      !! New solution vector.
   real(rk), intent(out) :: ydotnew(neq)
      !! New derivative of solution vector.
   real(rk), intent(inout) :: pwk(*)
      !! Real work array used by `psol`.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   real(rk), intent(in) :: rlx
      !! Factor for restricting update size in [[dcnstr]].
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   real(rk), parameter :: alpha = 1e-4_rk
   integer :: ierr, ivar, kprint
   real(rk) :: f1norm, f1normp, fnormp, ratio, ratio1, rl, rlmin, slpi, tau
   character(len=80) :: msg

   kprint = iwm(iwloc_kprint)
   f1norm = (fnorm**2)/2
   ratio = one

   if (kprint >= 2) then
      msg = '------ IN ROUTINE DLINSK-- PNRM = (R1)'
      call xerrwd(msg, 38, 921, 0, 0, 0, 0, 1, pnorm, zero)
   end if
   tau = pnorm
   rl = one

   ! Check for violations of the constraints, if any are imposed.
   ! If any violations are found, the step vector P is rescaled, and the
   ! constraint check is repeated, until no violations are found.
   if (icnflg /= 0) then
      do
         call dyypnw(neq, y, ydot, cj, rl, p, icopt, idy, ynew, ydotnew)
         call dcnstr(neq, y, ynew, icnstr, tau, rlx, iret, ivar)
         if (iret /= 1) exit
         ratio1 = tau/pnorm
         ratio = ratio*ratio1
         p = p*ratio1
         pnorm = tau
         if (kprint >= 2) then
            msg = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)'
            call xerrwd(msg, 50, 922, 0, 1, ivar, 0, 1, pnorm, zero)
         end if
         if (pnorm <= stptol) then
            iret = 1
            return
         end if
      end do
   end if

   slpi = -2*f1norm*ratio
   rlmin = stptol/pnorm
   if ((lsoff == 0) .and. (kprint >= 2)) then
      msg = '------ MIN. LAMBDA = (R1)'
      call xerrwd(msg, 25, 923, 0, 0, 0, 0, 1, rlmin, zero)
   end if

   ! Begin iteration to find RL value satisfying alpha-condition.
   ! Update YNEW and YDOTNEW, then compute norm of new scaled residual and
   ! perform alpha condition test.
   do
      call dyypnw(neq, y, ydot, cj, rl, p, icopt, idy, ynew, ydotnew)
      call dfnrmk(neq, ynew, t, ydotnew, savr, r, cj, tscale, wght, &
                  sqrtn, rsqrtn, res, ires, psol, 0, ierr, fnormp, epslin, &
                  rwp, iwp, pwk, rpar, ipar)
      iwm(iwloc_nre) = iwm(iwloc_nre) + 1
      if (ires >= 0) iwm(iwloc_nps) = iwm(iwloc_nps) + 1
      if (ires /= 0 .or. ierr /= 0) then
         iret = 2
         return
      end if

      if (lsoff == 1) goto 150

      f1normp = (fnormp**2)/2
      if (kprint >= 2) then
         msg = '------ LAMBDA = (R1)'
         call xerrwd(msg, 20, 924, 0, 0, 0, 0, 1, rl, zero)
         msg = '------ NORM(F1) = (R1),  NORM(F1NEW) = (R2)'
         call xerrwd(msg, 43, 925, 0, 0, 0, 0, 2, f1norm, f1normp)
      end if

      if (f1normp > f1norm + alpha*slpi*rl) goto 200

      ! Alpha-condition is satisfied, or linesearch is turned off.
      ! Copy YNEW,YDOTNEW to Y,YPRIME and return.
150   continue
      iret = 0
      y = ynew
      ydot = ydotnew
      fnorm = fnormp
      if (kprint >= 1) then
         msg = '------ LEAVING ROUTINE DLINSK, FNRM = (R1)'
         call xerrwd(msg, 42, 926, 0, 0, 0, 0, 1, fnorm, zero)
      end if
      return

      ! Alpha-condition not satisfied. Perform backtrack to compute new RL
      ! value. If RL is less than RLMIN, i.e. no satisfactory YNEW,YDOTNEW can
      ! be found sufficiently distinct from Y,YDOT, then return IRET = 1.
200   continue
      if (rl < rlmin) then
         iret = 1
         return
      end if

      rl = rl/2
   end do

end subroutine dlinsk

subroutine dfnrmk( &
   neq, y, t, ydot, savr, r, cj, tscale, wght, &
   sqrtn, rsqrtn, res, ires, psol, irin, ierr, &
   rnorm, epslin, rwp, iwp, wk, rpar, ipar)
!! This routine calculates the scaled preconditioned norm of the nonlinear function used
!! in the nonlinear iteration for obtaining consistent initial conditions. Specifically,
!! it calculates the weighted root-mean-square norm of the vector:
!!
!! $$  r = P^{-1} G(t,y,\dot{y}) $$
!!
!! where \(P\) is the preconditioner matrix and \(G\) is the DAE equation vector.

   use daskr_kinds, only: rk, zero
   use blas_interfaces, only: dcopy, dscal
   use daskr
   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(in) :: y(neq)
      !! Solution vector.
   real(rk), intent(in) :: t
      !! Independent variable.
   real(rk), intent(in) :: ydot(neq)
      !! Derivative of solution vector after successful step.
   real(rk), intent(inout) :: savr(neq)
      !! Saved residual vector.
   real(rk), intent(out) :: r(neq)
      !! Result vector.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: tscale ! @todo: what is "t?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: wght(neq)
      !! Scaling factors.
   real(rk), intent(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   integer, intent(in) :: irin
      !! Flag indicating whether the current residual vector is input in `savr`.
      !! `0`: it is not.
      !! `1`: it is.
   integer, intent(out) :: ierr
      !! Error flag from `psol`.
   real(rk), intent(out) :: rnorm
      !! Weighted root-mean-square norm of `r`.
   real(rk), intent(in) :: epslin
      !! Tolerance for linear system.
   real(rk), intent(inout) :: rwp(*)
      !! Real work array used by preconditioner `psol`.
   integer, intent(inout) :: iwp(*)
      !! Integer work array used by preconditioner `psol`.
   real(rk), intent(inout) :: wk(neq)
      !! Real work array used by `psol`.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   real(rk), external :: ddwnrm ! @todo: remove this once inside module

   ! Call RES routine if IRIN = 0.
   if (irin == 0) then
      ires = 0
      call res(t, y, ydot, cj, savr, ires, rpar, ipar)
      ! @todo: count missing here, or is it counted in the caller?
      if (ires < 0) return
   end if

   ! Apply inverse of left preconditioner to vector R.
   ! First scale WT array by 1/sqrt(N), and undo scaling afterward.
   call dcopy(neq, savr, 1, r, 1)
   call dscal(neq, rsqrtn, wght, 1)
   ierr = 0
   call psol(neq, t, y, ydot, savr, wk, cj, wght, rwp, iwp, r, epslin, ierr, rpar, ipar)
   ! @todo: count missing here, or is it counted in the caller?
   call dscal(neq, sqrtn, wght, 1)
   if (ierr /= 0) return

   ! Calculate norm of R.
   rnorm = ddwnrm(neq, r, wght, rpar, ipar)
   if (tscale > zero) rnorm = rnorm*tscale*abs(cj)

end subroutine dfnrmk

subroutine dnedk( &
   t, y, ydot, neq, res, jack, psol, &
   h, wt, jstart, idid, rpar, ipar, phi, gama, savr, delta, e, &
   rwm, iwm, cj, cjold, cjlast, s, uround, epli, sqrtn, rsqrtn, &
   epscon, jcalc, jflag, 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 matrix-free 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 after successful step.
   integer, intent(in) :: neq
      !! Problem size.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(jack_t) :: jack
      !! User-defined Jacobian routine.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   real(rk), intent(in) :: h
      !! Step size.
   real(rk), intent(inout) :: wt(neq)
      !! 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(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(out) :: cjold
      !! Value of `cj` as of the last call to [[ditmd]]. Accounts for changes in `cj`
      !! needed to decide whether to call [[ditmd]].
   real(rk), intent(in) :: cjlast
      !! Previous value of `cj`.
   real(rk), intent(out) :: s
      !! Convergence factor for the Newton iteration. See [[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 (not used).
   real(rk), intent(in) :: epli ! @note: not the same as 'epslin'!
      !! Convergence test constant for the iterative solution of linear systems. See
      !! [[daskr]] for details.
   real(rk), intent(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   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 [[ditmd]] to update the Jacobian matrix.
      !!  `0`: Jacobian matrix is up-to-date.
      !!  `1`: Jacobian matrix is outdated, but [[ditmd]] will not be called unless `jcalc` is set to `-1`.
   integer, intent(in) :: jflag
      !! Flag indicating whether a `jac` routine is supplied by the user.
   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.
      !! `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 = 0, maxit = 4
   real(rk), parameter :: xrate = 0.25_rk

   integer :: iernew, ierpj, iersl, iertyp, ires, j, liwp, lrwp
   real(rk) :: deltanorm, epslin, temp1, temp2, tolnew
   real(rk), external :: ddwnrm ! @todo: remove this once inside module

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

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

   ! Perform all other initializations.
   iernls = 0
   lrwp = iwm(iwloc_lrwp)
   liwp = iwm(iwloc_liwp)

   ! Decide whether to update the preconditioner.
   if (jflag /= 0) then
      temp1 = (one - xrate)/(one + xrate)
      temp2 = 1/temp1
      if ((cj/cjold < temp1) .or. (cj/cjold > temp2)) jcalc = -1
      if (cj /= cjlast) s = 1e2_rk
   else
      jcalc = 0
   end if

   ! Looping point for updating preconditioner with current stepsize.
300 continue

   ! Initialize all error flags to zero.
   ierpj = 0
   ires = 0
   iersl = 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

330 continue
   epslin = epli*epscon
   tolnew = epslin

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

   ! If indicated, update the preconditioner.
   ! Set JCALC to 0 as an indicator that this has been done.
   if (jcalc == -1) then
      jcalc = 0
      call jack(res, ires, neq, t, y, ydot, wt, delta, e, h, cj, rwm(lrwp), iwm(liwp), ierpj, rpar, ipar)
      iwm(iwloc_nje) = iwm(iwloc_nje) + 1
      cjold = cj
      s = 1e2_rk
      if (ires < 0) goto 380
      if (ierpj /= 0) goto 380
   end if

   ! Call the nonlinear Newton solver.
   call dnsk(t, y, ydot, neq, res, psol, wt, rpar, ipar, savr, &
             delta, e, rwm, iwm, cj, sqrtn, rsqrtn, epslin, epscon, &
             s, temp1, tolnew, muldel, maxit, ires, iersl, iernew)

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

   if (iernew /= 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 == 0) goto 390
   delta = min(y, zero)
   deltanorm = ddwnrm(neq, delta, wt, rpar, ipar)
   if (deltanorm > epscon) goto 380
   e = e - delta
   goto 390

   ! Exits from nonlinear solver.
   ! No convergence with current preconditioner.
   ! Compute IERNLS and IDID accordingly.
380 continue
   if ((ires <= -2) .or. (iersl < 0) .or. (iertyp /= 0)) then
      iernls = -1
      if (ires <= -2) idid = -11
      if (iersl < 0) idid = -13
      if (iertyp /= 0) idid = -15
   else
      iernls = 1
      if (ires == -1) idid = -10
      if (ierpj /= 0) idid = -5
      if (iersl > 0) idid = -14
   end if

390 continue
   jcalc = 1

end subroutine dnedk

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

subroutine dslvk( &
   neq, y, t, ydot, savr, x, ewt, rwm, iwm, res, ires, psol, &
   iersl, cj, epslin, sqrtn, rsqrtn, rhok, rpar, ipar)
!! This routine uses a restart algorithm and interfaces to [[dspigm]] for the solution of
!! the linear system arising from a Newton iteration.

   use daskr_kinds, only: rk, zero
   use daskr
   use blas_interfaces, only: dscal, dcopy
   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(in) :: y(neq)
      !! Current dependent variables.
   real(rk), intent(in) :: t
      !! Current independent variable.
   real(rk), intent(in) :: ydot(neq)
      !! Current derivatives of dependent variables.
   real(rk), intent(in) :: savr(neq)
      !! Current residual evaluated at `(t, y, ydot)`.
   real(rk), intent(inout) :: x(neq)
      !! On entry, the right-hand side vector of the linear system to be solved,
      !! and, on exit, the solution vector.
   real(rk), intent(inout) :: ewt(neq)
      !! Nonzero elements of the diagonal scaling matrix.
   real(rk), intent(inout) :: rwm(*)
      !! Real work space containing data for the algorithm (Krylov basis vectors,
      !! Hessenberg matrix, etc.).
   integer, intent(inout) :: iwm(*)
      !! Integer work space containing data for the algorithm.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   integer, intent(out) :: iersl
      !! Error flag.
      !!  `0`: no trouble occurred (or user `res` routine returned `ires < 0`).
      !!  `1`: iterative method failed to converge ([[dspigm]] returned `iflag > 0`).
      !! `-1`: nonrecoverable error in the iterative solver, and an error exit will occur.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: epslin
      !! Tolerance for linear system solver.
   real(rk), intent(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   real(rk), intent(out) :: rhok
      !! Weighted norm of the final preconditioned residual.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   integer, save :: irst = 1 ! @todo: move this to iwork, to get rid of save
   integer :: i, iflag, kmp, ldl, lhes, lgmr, liwp, lq, lr, lv, lwk, lrwp, lz, &
              maxl, maxlp1, miter, ncfl, nli, nps, npsol, nre, nres, nrmax, nrsts

   liwp = iwm(iwloc_liwp)
   nli = iwm(iwloc_nli)
   nps = iwm(iwloc_nps)
   ncfl = iwm(iwloc_ncfl)
   nre = iwm(iwloc_nre)
   lrwp = iwm(iwloc_lrwp)
   maxl = iwm(iwloc_maxl)
   kmp = iwm(iwloc_kmp)
   nrmax = iwm(iwloc_nrmax)
   miter = iwm(iwloc_miter)
   iersl = 0
   ires = 0

   ! Use a restarting strategy to solve the linear system P*X = -F. Parse the work vector,
   ! and perform initializations. Note that zero is the initial guess for X.
   maxlp1 = maxl + 1
   lv = 1
   lr = lv + neq*maxl
   lhes = lr + neq + 1
   lq = lhes + maxl*maxlp1
   lwk = lq + 2*maxl
   ldl = lwk + min(1, maxl - kmp)*neq
   lz = ldl + neq
   call dscal(neq, rsqrtn, ewt, 1)
   call dcopy(neq, x, 1, rwm(lr), 1)
   x = zero

   ! Top of loop for the restart algorithm. Initial pass approximates X and sets up a
   ! transformed system to perform subsequent restarts to update X. NRSTS is initialized
   ! to -1, because restarting does not occur until after the first pass.
   ! Update NRSTS; conditionally copy DL to R; call the DSPIGM algorithm to solve A*Z = R;
   ! updated counters; update X with the residual solution.
   ! Note: if convergence is not achieved after NRMAX restarts, then the linear solver is
   ! considered to have failed.
   iflag = 1
   nrsts = -1
   do while ((iflag == 1) .and. (nrsts < nrmax) .and. (ires == 0))
      nrsts = nrsts + 1
      if (nrsts > 0) call dcopy(neq, rwm(ldl), 1, rwm(lr), 1)
      call dspigm(neq, t, y, ydot, savr, rwm(lr), ewt, maxl, &
                  kmp, epslin, cj, res, ires, nres, psol, npsol, rwm(lz), rwm(lv), &
                  rwm(lhes), rwm(lq), lgmr, rwm(lrwp), iwm(liwp), rwm(lwk), &
                  rwm(ldl), rhok, iflag, irst, nrsts, rpar, ipar)
      nli = nli + lgmr
      nps = nps + npsol
      nre = nre + nres
      do i = 1, neq
         x(i) = x(i) + rwm(lz + i - 1)
      end do
   end do

   ! The restart scheme is finished. Test IRES and IFLAG to see if convergence was not
   ! achieved, and set flags accordingly.
   if (ires < 0) then
      ncfl = ncfl + 1
   elseif (iflag /= 0) then
      ncfl = ncfl + 1
      if (iflag > 0) iersl = 1
      if (iflag < 0) iersl = -1
   end if

   ! Update IWM with counters, rescale EWT, and return.
   iwm(iwloc_nli) = nli
   iwm(iwloc_nps) = nps
   iwm(iwloc_ncfl) = ncfl
   iwm(iwloc_nre) = nre
   call dscal(neq, sqrtn, ewt, 1)

end subroutine dslvk

subroutine dspigm( &
   neq, t, y, ydot, savr, r, wght, maxl, &
   kmp, epslin, cj, res, ires, nres, psol, npsol, z, v, &
   hes, q, lgmr, rwp, iwp, wk, dl, rhok, iflag, irst, nrsts, &
   rpar, ipar)
!! This routine solves the linear system \(A z = r\) using a scaled preconditioned version
!! of the generalized minimum residual method. An initial guess of \(z=0\) is assumed.

   use daskr_kinds, only: rk, zero, one
   use daskr
   use blas_interfaces, only: daxpy, dcopy, dnrm2, dscal

   implicit none

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(in) :: t
      !! Current independent variable.
   real(rk), intent(in) :: y(neq)
      !! Current dependent variables.
   real(rk), intent(in) :: ydot(neq)
      !! Current derivatives of dependent variables.
   real(rk), intent(in) :: savr(neq)
      !! Current residual evaluated at `(t, y, ydot)`.
   real(rk), intent(inout) :: r(neq)
      !! On entry, the right hand side vector. Also used as workspace when computing
      !! the final approximation and will therefore be destroyed. `r` is the same as
      !! `v(:,maxl+1)` in the call to this routine.
   real(rk), intent(in) :: wght(neq)
      !! Nonzero elements of the diagonal scaling matrix.
   integer, intent(in) :: maxl
      !! Maximum allowable order of the matrix `hes`.
   integer, intent(in) :: kmp
      !! Number of previous vectors the new vector `vnew` must be made orthogonal to
      !! (`kmp <= maxl`).
   real(rk), intent(in) :: epslin
      !! Tolerance on residuals \(Az - r\) in weighted rms norm.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   integer, intent(inout) :: nres
      !! Number of calls to `res`.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   integer, intent(inout) :: npsol
      !! Number of calls to `psol`.
   real(rk), intent(out) :: z(neq)
      !! Final computed approximation to the solution of the system \(Az = r\).
   real(rk), intent(out) :: v(neq, *)
      !! Matrix of shape `(neq,lgmr+1)` containing the `lgmr` orthogonal vectors `v(:,1)`
      !! to `v(:,lgmr)`.
   real(rk), intent(out) :: hes(maxl + 1, maxl)
      !! Upper triangular factor of the QR decomposition of the upper Hessenberg matrix
      !! whose entries are the scaled inner-products of `A*v(:,i)` and `v(:,k)`.
   real(rk), intent(out) :: q(2*maxl)
      !! Components of the Givens rotations used in the QR decomposition of `hes`. It is
      !! loaded in [[dheqr]] and used in [[dhels]].
   integer, intent(out) :: lgmr
      !! The number of iterations performed and the current order of the upper
      !! Hessenberg matrix `hes`.
   real(rk), intent(inout) :: rwp(*)
      !! Real work array used by preconditioner `psol`.
   integer, intent(inout) :: iwp(*)
      !! Integer work array used by preconditioner `psol`.
   real(rk), intent(inout) :: wk(*)
      !! Real work array used by [[datv]] and `psol`.
   real(rk), intent(inout) :: dl(*)
      !! Real work array used for calculation of the residual norm `rhok` when the
      !! method is incomplete (`kmp < maxl`) and/or when using restarting.
   real(rk), intent(out) :: rhok
      !! Weighted norm of the final preconditioned residual.
   integer, intent(out) :: iflag
      !! Error flag.
      !! `0`: convergence in `lgmr` iterations, `lgmr <= maxl`.
      !! `1`: convergence test did not pass in `maxl` iterations, but the new
      !! residual norm (`rho`) is less than the old residual norm (`rnrm`), and so `z`
      !! is computed.
      !! `2`: convergence test did not pass in `maxl` iterations, new residual
      !! norm (`rho`) is greater than or equal to old residual norm (`rnrm`), and the
      !! initial guess, `z = 0`, is returned.
      !! `3`: recoverable error in `psol` caused by the preconditioner being out of date.
      !! `-1`: unrecoverable error in `psol`.
   integer, intent(in) :: irst
      !! Restarting flag. If `irst > 0`, then restarting is being performed.
   integer, intent(in) :: nrsts
      !! Number of restarts on the current call to [[dspigm]]. If `nrsts > 0`, then the
      !! residual `r` is already scaled, and so scaling of `r` is not necessary.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   integer :: i, ierr, info, ip1, i2, k, ll, llp1, maxlm1, maxlp1
   real(rk) :: c, dlnorm, prod, rho, rnorm, s, snormw, tem

   ierr = 0
   iflag = 0
   lgmr = 0
   npsol = 0
   nres = 0

   ! The initial guess for Z is 0. The initial residual is therefore
   ! the vector R. Initialize Z to 0.
   z = zero

   ! Apply inverse of left preconditioner to vector R if NRSTS == 0.
   ! Form V(:,1), the scaled preconditioned right hand side.
   if (nrsts == 0) then
      call psol(neq, t, y, ydot, savr, wk, cj, wght, rwp, iwp, r, epslin, ierr, rpar, ipar)
      npsol = 1
      if (ierr /= 0) goto 300
      v(:, 1) = r*wght
   else
      v(:, 1) = r
   end if

   ! Calculate norm of scaled vector V(:,1) and normalize it
   ! If, however, the norm of V(:,1) (i.e. the norm of the preconditioned
   ! residual) is <= EPSLIN, then return with Z=0.
   rnorm = dnrm2(neq, v, 1)
   if (rnorm <= epslin) then
      rhok = rnorm
      return
   end if
   tem = one/rnorm
   call dscal(neq, tem, v(1, 1), 1)

   ! Zero out the HES array.
   hes = zero

   ! Main loop to compute the vectors V(:,2) to V(:,MAXL).
   ! The running product PROD is needed for the convergence test.
   prod = one
   maxlp1 = maxl + 1
   do ll = 1, maxl
      lgmr = ll

      ! Call routine DATV to compute VNEW = ABAR*V(LL), where ABAR is
      ! the matrix A with scaling and inverse preconditioner factors applied.
      call datv(neq, y, t, ydot, savr, v(1, ll), wght, z, &
                res, ires, psol, v(1, ll + 1), wk, rwp, iwp, cj, epslin, &
                ierr, nres, npsol, rpar, ipar)
      if (ires < 0) return
      if (ierr /= 0) goto 300

      ! Call routine DORTH to orthogonalize the new vector VNEW = V(:,LL+1).
      call dorth(v(1, ll + 1), v, hes, neq, ll, maxlp1, kmp, snormw)
      hes(ll + 1, ll) = snormw

      ! Call routine DHEQR to update the factors of HES.
      call dheqr(hes, maxlp1, ll, q, info, ll)
      if (info == ll) goto 120

      ! Update RHO, the estimate of the norm of the residual R - A*ZL.
      ! If KMP < MAXL, then the vectors V(:,1),...,V(:,LL+1) are not
      ! necessarily orthogonal for LL > KMP.  The vector DL must then
      ! be computed, and its norm used in the calculation of RHO.
      prod = prod*q(2*ll)
      rho = abs(prod*rnorm)
      if ((ll > kmp) .and. (kmp < maxl)) then
         if (ll == kmp + 1) then
            call dcopy(neq, v(1, 1), 1, dl, 1)
            do i = 1, kmp
               ip1 = i + 1
               i2 = i*2
               s = q(i2)
               c = q(i2 - 1)
               do k = 1, neq
                  dl(k) = s*dl(k) + c*v(k, ip1)
               end do
            end do
         end if
         s = q(2*ll)
         c = q(2*ll - 1)/snormw
         llp1 = ll + 1
         do k = 1, neq
            dl(k) = s*dl(k) + c*v(k, llp1)
         end do
         dlnorm = dnrm2(neq, dl, 1)
         rho = rho*dlnorm
      end if

      ! Test for convergence. If passed, compute approximation ZL.
      ! If failed and LL < MAXL, then continue iterating.
      if (rho <= epslin) goto 200
      if (ll == maxl) goto 100

      ! Rescale so that the norm of V(1,LL+1) is one.
      tem = one/snormw
      call dscal(neq, tem, v(1, ll + 1), 1)

   end do

100 continue
   if (rho < rnorm) goto 150

120 continue
   iflag = 2
   z = zero
   return

150 continue
   ! The tolerance was not met, but the residual norm was reduced.
   ! If performing restarting (IRST > 0) calculate the residual vector
   ! RL and store it in the DL array.  If the incomplete version is
   ! being used (KMP < MAXL) then DL has already been calculated.
   iflag = 1
   if (irst > 0) then

      if (kmp == maxl) then
         !  Calculate DL from the V(I)'s.
         call dcopy(neq, v(1, 1), 1, dl, 1)
         maxlm1 = maxl - 1
         do i = 1, maxlm1
            ip1 = i + 1
            i2 = i*2
            s = q(i2)
            c = q(i2 - 1)
            do k = 1, neq
               dl(k) = s*dl(k) + c*v(k, ip1)
            end do
         end do
         s = q(2*maxl)
         c = q(2*maxl - 1)/snormw
         do k = 1, neq
            dl(k) = s*dl(k) + c*v(k, maxlp1)
         end do
      end if

      ! Scale DL by RNRM*PROD to obtain the residual RL.
      tem = rnorm*prod
      call dscal(neq, tem, dl, 1)
   end if

   ! Compute the approximation ZL to the solution.
   ! Since the vector Z was used as work space, and the initial guess
   ! of the Newton correction is zero, Z must be reset to zero.
200 continue
   ll = lgmr
   llp1 = ll + 1
   r(1:llp1) = zero
   r(1) = rnorm
   call dhels(hes, maxlp1, ll, q, r)
   z = zero
   do i = 1, ll
      call daxpy(neq, r(i), v(1, i), 1, z, 1)
   end do
   z = z/wght
   ! Load RHO into RHOK.
   rhok = rho
   return

   ! This block handles error returns forced by routine PSOL.
300 continue
   if (ierr < 0) iflag = -1
   if (ierr > 0) iflag = 3

end subroutine dspigm

subroutine datv( &
   neq, y, t, ydot, savr, v, wght, ydottemp, res, &
   ires, psol, z, vtemp, rwp, iwp, cj, epslin, ierr, nres, npsol, rpar, ipar)
!! This routine computes the matrix-vector product:
!!
!! $$ z = D^{-1} P^{-1} (\partial F / \partial y) D v $$
!!
!! where \(F = G(t, y, c_J(y - a))\), \(c_J\) is a scalar proportional to \(1/h\), and \(a\)
!! involves the past history of \(y\). The quantity \(c_J(y - a)\) is an approximation to
!! the first derivative of \(y\) and is stored in \(\dot{y}\).
!!
!! \(D\) is a diagonal scaling matrix, and \(P\) is the left preconditioning matrix. \(v\)
!! is assumed to have L-2 norm equal to 1. The product is stored in \(z\) and is computed by
!! means of a difference quotient, a call to `res`, and one call to `psol`.

   use daskr_kinds, only: rk
   use daskr

   integer, intent(in) :: neq
      !! Problem size.
   real(rk), intent(in) :: y(neq)
      !! Current dependent variables.
   real(rk), intent(in) :: t
      !! Current independent variable.
   real(rk), intent(in) :: ydot(neq)
      !! Current derivatives of dependent variables.
   real(rk), intent(in) :: savr(neq)
      !! Current residual evaluated at `(t, y, ydot)`.
   real(rk), intent(in) :: v(neq)
      !! Orthonormal vector (can be the same as `z`).
   real(rk), intent(in) :: wght(neq)
      !! Scaling factors. `1/wght(i)` are the diagonal elements of the matrix \(D\).
   real(rk), intent(out) :: ydottemp(neq)
      !! Work array used to store the incremented value of `ydot`.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   integer, intent(out) :: ires
      !! Error flag from `res`.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   real(rk), intent(out) :: z(neq)
      !! Desired scaled matrix-vector product.
   real(rk), intent(out) :: vtemp(neq)
      !! Work array used to store the unscaled version of `v`.
   real(rk), intent(inout) :: rwp(*)
      !! Real work array used by preconditioner `psol`.
   integer, intent(inout) :: iwp(*)
      !! Integer work array used by preconditioner `psol`.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(in) :: epslin
      !! Tolerance for linear system.
   integer, intent(out) :: ierr
      !! Error flag from `psol`.
   integer, intent(inout) :: nres
      !! Number of calls to `res`.
   integer, intent(inout) :: npsol
      !! Number of calls to `psol`.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.

   ! Set VTEMP = D * V.
   vtemp = v/wght

   ! Store Y in Z and increment Z by VTEMP.
   ! Store YDOT in YDOTTEMP and increment YDOTTEMP by VTEM*CJ.
   ydottemp = ydot + vtemp*cj
   z = y + vtemp

   ! Call RES with incremented Y, YDOT arguments
   ! stored in Z, YDOTTEMP. VTEMP is overwritten with new residual.
   ires = 0
   call res(t, z, ydottemp, cj, vtemp, ires, rpar, ipar)
   nres = nres + 1
   if (ires < 0) return

   ! Set Z = (dF/dY) * VBAR using difference quotient.
   ! (VBAR is old value of VTEMP before calling RES)
   z = vtemp - savr

   ! Apply inverse of left preconditioner to Z.
   ierr = 0
   call psol(neq, t, y, ydot, savr, ydottemp, cj, wght, rwp, iwp, z, epslin, ierr, rpar, ipar)
   npsol = npsol + 1
   if (ierr /= 0) return

   ! Apply D-inverse to Z
   z = z*wght

end subroutine datv

pure subroutine dorth(vnew, v, hes, n, ll, ldhes, kmp, snormw)
!! This routine orthogonalizes the vector `vnew` against the previous `kmp` vectors in the
!! `v` matrix. It uses a modified Gram-Schmidt orthogonalization procedure with conditional
!! reorthogonalization.

   use daskr_kinds, only: rk, zero
   use blas_interfaces, only: daxpy, ddot, dnrm2
   implicit none

   real(rk), intent(inout) :: vnew(n)
      !! On entry, vector containing a scaled product of the Jacobian and the vector `v(:,ll)`.
      !! On return, the new vector orthogonal to `v(:,i0)`, where `i0 = max(1, ll - kmp + 1)`.
   real(rk), intent(in) :: v(n, ll)
      !! Matrix containing the previous `ll` orthogonal vectors `v(:,1)` to `v(:,ll)`.
   real(rk), intent(inout) :: hes(ldhes, ll)
      !! On entry, an upper Hessenberg matrix of shape `(ll, ll)` containing in `hes(i,k)`,
      !! for `k < ll`, the scaled inner products of `a*v(:,k)` and `v(:,i)`.
      !! On return, an upper Hessenberg matrix with column `ll` filled in with the scaled
      !! inner products of `a*v(:,ll)` and `v(:,i)`.
   integer, intent(in) :: n
      !! Order of the matrix `a`.
   integer, intent(in) :: ll
      !! Current order of the matrix `hes`.
   integer, intent(in) :: ldhes
      !! Leading dimension of `hes`.
   integer, intent(in) :: kmp
      !! Number of previous vectors the new vector `vnew` must be made orthogonal to
      !! (`kmp <= maxl`).
   real(rk), intent(out) :: snormw
      !! L-2 norm of `vnew`.

   integer :: i, i0
   real(rk) :: arg, sumdsq, tem, vnorm

   ! Get norm of unaltered VNEW for later use.
   vnorm = dnrm2(n, vnew, 1)

   ! Do Modified Gram-Schmidt on VNEW = A*V(LL).
   ! Scaled inner products give new column of HES.
   ! Projections of earlier vectors are subtracted from VNEW.
   i0 = max(1, ll - kmp + 1)
   do i = i0, ll
      hes(i, ll) = ddot(n, v(1, i), 1, vnew, 1)
      tem = -hes(i, ll)
      call daxpy(n, tem, v(1, i), 1, vnew, 1)
   end do

   ! Compute SNORMW = norm of VNEW.
   ! If VNEW is small compared to its input value (in norm), then
   ! Reorthogonalize VNEW to V(*,1) through V(*,LL).
   ! Correct if relative correction exceeds 1000*(unit roundoff).
   ! Finally, correct SNORMW using the dot products involved.
   snormw = dnrm2(n, vnew, 1)
   if (vnorm + snormw/1000 /= vnorm) return ! @todo: fix this comparison

   sumdsq = zero
   do i = i0, ll
      tem = -ddot(n, v(1, i), 1, vnew, 1)
      if (hes(i, ll) + tem/1000 == hes(i, ll)) cycle ! @todo: fix this comparison
      hes(i, ll) = hes(i, ll) - tem
      call daxpy(n, tem, v(1, i), 1, vnew, 1)
      sumdsq = sumdsq + tem**2
   end do
   if (sumdsq == zero) return

   arg = max(zero, snormw**2 - sumdsq)
   snormw = sqrt(arg)

end subroutine dorth

pure subroutine dheqr(a, lda, n, q, info, ijob)
!! This routine performs a QR decomposition of an upper Hessenberg matrix `a` using Givens
!! rotations. There are two options available:
!!
!! 1. Performing a fresh decomposition;
!! 2. Updating the QR factors by adding a row and a column to the matrix `a`.

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

   real(rk), intent(inout) :: a(lda, n)
      !! On entry, the Hessenberg matrix to be decomposed. On return, the upper triangular
      !! matrix R from the QR decomposition of `a`.
   integer, intent(in) :: lda
      !! Leading dimension of `a`.
   integer, intent(in) :: n
      !! Number of columns in `a` (originally a matrix with shape `(n+1, n)`).
   real(rk), intent(out) :: q(2*n)
      !! Coefficients of the Givens rotations used in decomposing `a`.
   integer, intent(out) :: info
      !! Info flag.
      !! `info = 0`, if successful.
      !! `info = k`, if `a(k,k) = 0`; this is not an error condition for this subroutine, but
      !! it does indicate that [[dhels]] will divide by zero if called.
   integer, intent(in) :: ijob
      !! Job flag.
      !! `ijob = 1`, means that a fresh decomposition of the matrix `a` is desired.
      !! `ijob >= 2`, means that the current decomposition of `a` will be updated by the
      !! addition of a row and a column.

   integer :: i, iq, j, k, km1, kp1, nm1
   real(rk) :: c, s, t, t1, t2

   ! A new factorization is desired.
   if (ijob == 1) then

      ! QR decomposition without pivoting.
      info = 0
      do k = 1, n
        
         km1 = k - 1
         kp1 = k + 1

         ! Compute Kth column of R.
         ! First, multiply the Kth column of A by the previous K-1 Givens rotations.
         if (km1 < 1) goto 20
         do j = 1, km1
            i = 2*(j - 1) + 1
            t1 = a(j, k)
            t2 = a(j + 1, k)
            c = q(i)
            s = q(i + 1)
            a(j, k) = c*t1 - s*t2
            a(j + 1, k) = s*t1 + c*t2
         end do

         ! Compute Givens components C and S.
20       continue
         iq = 2*km1 + 1
         t1 = a(k, k)
         t2 = a(kp1, k)
         if (t2 /= zero) goto 30
         c = one
         s = zero
         goto 50

30       continue
         if (abs(t2) < abs(t1)) goto 40
         t = t1/t2
         s = -one/sqrt(one + t*t)
         c = -s*t
         goto 50

40       continue
         t = t2/t1
         c = one/sqrt(one + t*t)
         s = -c*t

50       continue
         q(iq) = c
         q(iq + 1) = s
         a(k, k) = c*t1 - s*t2
         if (a(k, k) == zero) info = k
      end do

      ! The old factorization of A will be updated.  A row and a column
      ! has been added to the matrix A.
      ! N by N-1 is now the old size of the matrix.
   else if (ijob >= 2) then

      nm1 = n - 1

      ! Multiply the new column by the N previous Givens rotations.
      do k = 1, nm1
         i = 2*(k - 1) + 1
         t1 = a(k, n)
         t2 = a(k + 1, n)
         c = q(i)
         s = q(i + 1)
         a(k, n) = c*t1 - s*t2
         a(k + 1, n) = s*t1 + c*t2
      end do

      ! Complete update of decomposition by forming last Givens rotation,
      ! and multiplying it times the column vector (A(N,N),A(NP1,N)).
      info = 0
      t1 = a(n, n)
      t2 = a(n + 1, n)
      if (t2 /= zero) goto 110
      c = one
      s = zero
      goto 130

110   continue
      if (abs(t2) < abs(t1)) goto 120
      t = t1/t2
      s = -one/sqrt(one + t*t)
      c = -s*t
      goto 130

120   continue
      t = t2/t1
      c = one/sqrt(one + t*t)
      s = -c*t

130   continue
      iq = 2*n - 1
      q(iq) = c
      q(iq + 1) = s
      a(n, n) = c*t1 - s*t2
      if (a(n, n) == zero) info = n

   end if

end subroutine dheqr

pure subroutine dhels(a, lda, n, q, b)
!! This routine solves the least squares problem:
!!
!!  $$ \min{\| b - A x \|^2} $$
!!
!! using the factors computed by [[dheqr]]. This is similar to the LINPACK routine [[dgesl]]
!! except that \(A\) is an upper Hessenberg matrix.

   use daskr_kinds, only: rk
   use blas_interfaces, only: daxpy
   implicit none

   real(rk), intent(in) :: a(lda, n)
      !! Output from [[dheqr]] which contains the upper triangular factor R in the QR
      !! decomposition of `a`.
   integer, intent(in) :: lda
      !! Leading dimension of `a`.
   integer, intent(in) :: n
      !! Number of columns in `a` (originally a matrix with shape `(n+1, n)`).
   real(rk), intent(in) :: q(2*n)
      !! Coefficients of the Givens rotations used in decomposing `a`.
   real(rk), intent(inout) :: b(n + 1)
      !! On entry, the right hand side vector. On return, the solution vector \(x\).

   integer :: iq, k, kb, kp1
   real(rk) :: c, s, t, t1, t2

   ! Minimize (B-A*X,B-A*X).
   ! First form Q*B.
   do k = 1, n
      kp1 = k + 1
      iq = 2*(k - 1) + 1
      c = q(iq)
      s = q(iq + 1)
      t1 = b(k)
      t2 = b(kp1)
      b(k) = c*t1 - s*t2
      b(kp1) = s*t1 + c*t2
   end do

   ! Now solve R*X = Q*B.
   do kb = 1, n
      k = n + 1 - kb
      b(k) = b(k)/a(k, k)
      t = -b(k)
      call daxpy(k - 1, t, a(1, k), 1, b(1), 1)
   end do

end subroutine dhels