banja Subroutine

   subroutine banja(res, ires, neq, t, y, yprime, rewt, savr, wk, h, cj, wp, iwp, ier, rpar, ipar)
   !! This subroutine generates a banded preconditioner matrix \(P\) that approximates the
   !! iteration matrix \(J = dG/dy + c_j dG/dy'\), where the DAE system is \(G(t,y,y') = 0\). 
   !! The band matrix \(P\) has half-bandwidths \(\mathrm{ml}\) and \(\mathrm{mu}\). It 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]]. 
   !! [[banja]] calls the LINPACK routine [[DGBFA]] to do an LU factorization of this matrix.
      external :: res
      integer, intent(inout) :: ires
        !! Output flag set by `res`. See `res` description in [[daskr]].
      integer, intent(in) :: neq
        !! Problem size.
      real(rk), intent(in) :: t
        !! Independent variable.
      real(rk), intent(inout) :: y(*)
        !! Current dependent variables.
      real(rk), intent(inout) :: yprime(*)
        !! Current derivatives of dependent variables.
      real(rk), intent(in) :: rewt(*)
        !! Vector of reciprocal error weights, used here for computing increments.
      real(rk), intent(in) :: savr(*)
        !! Current residual evaluated at `(t, y, yprime)`.
      real(rk), intent(in) :: wk(*)
        !! Real work space of length `neq`.
      real(rk), intent(in) :: h
        !! Current step size.
      real(rk), intent(in) :: cj
        !! Scalar proportional to `1/h`.
      real(rk), intent(inout) :: wp(*)
        !! 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) :: ier
        !! Output flag: `ier > 0` if P is singular, and `ier = 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.

      external :: dgbfa, d1mach

      real(rk) :: del, delinv, squr, uround, d1mach
      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.
      uround = d1mach(4)
      squr = sqrt(uround)

      ! 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 YPRIME elements.
      lenp = (2*ml + mu + 1)*neq
      msave = (neq/mband) + 1
      isave = lenp
      ipsave = isave + msave

      ! Initialize error flags.
      ier = 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
            wp(isave + k) = y(n)
            wp(ipsave + k) = yprime(n)
            del = squr*max(abs(y(n)), abs(h*yprime(n)), abs(one/rewt(n)))
            del = sign(del, h*yprime(n))
            del = (y(n) + del) - y(n)
            y(n) = y(n) + del
            yprime(n) = yprime(n) + cj*del
         end do
         
         call res(t, y, yprime, cj, wk, ires, rpar, ipar)
         
         if (ires < 0) return
         
         do n = j, neq, mband
            k = (n - j)/mband + 1
            y(n) = wp(isave + k)
            yprime(n) = wp(ipsave + k)
            del = squr*max(abs(y(n)), abs(h*yprime(n)), abs(one/rewt(n)))
            del = sign(del, h*yprime(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
               wp(ii + i) = (wk(i) - savr(i))*delinv
            end do
         end do
      
      end do

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

   end subroutine banja