jac_banpre Subroutine

   subroutine jac_banpre( &
      res, ires, neq, t, y, ydot, rewt, savr, wk, h, cj, rwp, iwp, ierr, 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) :: ires
        !! 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
        !! Error flag: `ierr > 0` if \(P\) is singular, and `ierr = 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, squround
      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.
      squround = 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 = 0
      ires = 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 = squround*max(abs(y(n)), abs(h*ydot(n)), abs(one/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, ires, rpar, ipar)
         
         if (ires < 0) return
         
         do n = j, neq, mband
            k = (n - j)/mband + 1
            y(n) = rwp(isave + k)
            ydot(n) = rwp(ipsave + k)
            del = squround*max(abs(y(n)), abs(h*ydot(n)), abs(one/rewt(n)))
            del = sign(del, h*ydot(n))
            del = (y(n) + del) - y(n)
            delinv = one/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)

   end subroutine jac_banpre