dslvk Subroutine

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