pjac_banpre Subroutine

   subroutine pjac_banpre( &
      res, ierr_res, neq, t, y, ydot, rewt, savr, wk, h, cj, rwp, iwp, ierr_pjac, rpar, ipar)
   !! This routine generates a banded preconditioner matrix \(P\) that approximates the
   !! Jacobian matrix \(J\). The banded matrix \(P\) has half-bandwidths \(\mathrm{ml}\)
   !! and \(\mathrm{mu}\), and is computed by making \(\mathrm{ml} + \mathrm{mu} + 1\) calls
   !! to the user's `res` routine and forming difference quotients, exactly as in the banded
   !! direct method option of [[daskr]]. Afterwards, this matrix is LU factorized by [[dgbfa]]
   !! from LINPACK and the factors are stored in the work arrays `rwp` and `iwp`.

      use linpack, only: dgbfa

      procedure(res_t) :: res
        !! User-defined residuals routine.
      integer, intent(out) :: ierr_res
        !! Error flag set by `res`.
      integer, intent(in) :: neq
        !! Problem size.
      real(rk), intent(in) :: t
        !! Current independent variable.
      real(rk), intent(inout) :: y(*)
        !! Current dependent variables.
      real(rk), intent(inout) :: ydot(*)
        !! Current derivatives of dependent variables.
      real(rk), intent(in) :: rewt(*)
        !! Reciprocal error weights for scaling `y` and `ydot`.
      real(rk), intent(inout) :: savr(*)
        !! Current residual evaluated at `(t, y, ydot)`.
      real(rk), intent(inout) :: wk(*)
        !! Real work space available to this subroutine.
      real(rk), intent(in) :: h
        !! Current step size.
      real(rk), intent(in) :: cj
        !! Scalar used in forming the system Jacobian.
      real(rk), intent(inout) :: rwp(*)
        !! Real work array for \(P\), etc. On output, it contains the LU decomposition of
        !! the banded approximation \(P\).
      integer, intent(inout) :: iwp(*)
        !! Integer work space for matrix pivot information.
      integer, intent(inout) :: ierr_pjac
        !! Error flag, `> 0` if \(P\) is singular, and `0` otherwise.
      real(rk), intent(inout) :: rpar(*)
        !! Real array used for communication between the calling program and external user
        !! routines.
      integer, intent(inout) :: ipar(*)
        !! Integer array used for communication between the calling program and external
        !! user routines. `ipar(1)` and `ipar(2)` must contain `ml` and `mu`, respectively.

      real(rk) :: del, delinv, sqrt_uround
      integer :: i, i1, i2, ii, ipsave, isave, j, k, lenp, mba, mband, meb1, meband, ml, &
                 msave, mu, n

      ! Set band parameters.
      ml = ipar(1)
      mu = ipar(2)
      mband = ml + mu + 1
      mba = min(mband, neq)
      meband = mband + ml
      meb1 = meband - 1

      ! Set the machine unit roundoff UROUND and SQRT(UROUND), used to set increments in
      ! the difference quotient procedure.
      sqrt_uround = sqrt(epsilon(one))

      ! Set pointers into WP. LENP is the length of the segment for P.
      ! Following that are two segments of size (NEQ/MBAND), with offsets
      ! ISAVE and IPSAVE, for temporary storage of Y and YDOT elements.
      lenp = (2*ml + mu + 1)*neq
      msave = (neq/mband) + 1
      isave = lenp
      ipsave = isave + msave

      ! Initialize error flags.
      ierr_pjac = 0
      ierr_res = 0

      ! Generate the banded approximate iteration matrix P using difference quotients on
      ! the results of calls to RES.
      do j = 1, mba

         do n = j, neq, mband
            k = (n - j)/mband + 1
            rwp(isave + k) = y(n)
            rwp(ipsave + k) = ydot(n)
            del = sqrt_uround*max(abs(y(n)), abs(h*ydot(n)), abs(1/rewt(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, wk, ierr_res, rpar, ipar)

         if (ierr_res < 0) return

         do n = j, neq, mband
            k = (n - j)/mband + 1
            y(n) = rwp(isave + k)
            ydot(n) = rwp(ipsave + k)
            del = sqrt_uround*max(abs(y(n)), abs(h*ydot(n)), abs(1/rewt(n)))
            del = sign(del, h*ydot(n))
            del = (y(n) + del) - y(n)
            delinv = 1/del
            i1 = max(1, n - mu)
            i2 = min(neq, n + ml)
            ii = n*meb1 - ml
            do i = i1, i2
               rwp(ii + i) = (wk(i) - savr(i))*delinv
            end do
         end do

      end do

      ! Do LU decomposition of the band matrix P.
      call dgbfa(rwp, meband, neq, ml, mu, iwp, ierr_pjac)

   end subroutine pjac_banpre