jac_rbgpre Subroutine

   subroutine jac_rbgpre(t, u, r0, rblock, r1, rewt, cj, bd, ipbd, ierr)
   !! This routine generates and preprocesses a block-diagonal preconditioner matrix, based on
   !! the part of the Jacobian corresponding to the reaction terms \(R\) of the problem, using
   !! block-grouping. It generates a matrix of the form \(c_J I_d - \partial R/ \partial u\).
   !! It calls [[dgefa]] from LINPACK to do the LU decomposition of each diagonal block. The
   !! computation of the diagonal blocks uses the mesh information in the module variables. One
   !! block per group is computed. The Jacobian elements are generated by difference quotients. 
   !! This routine calls a user routine `rblock` to compute the block `(jx,jy)` of \(R\).
      
      use linpack, only: dgefa

      real(rk), intent(in) :: t
         !! Current independent variable.
      real(rk), intent(inout) :: u(*)
         !! Current dependent variables.
      real(rk), intent(in) :: r0(*)
         !! Current value of the vector \(R(t,u)\).
      interface
         subroutine rblock(t, jx, jy, uxy, rxy)
         !! User routine that computes a single block of \(R\).
            import :: rk
            real(rk), intent(in) :: t
               !! Current independent variable.
            integer, intent(in) :: jx
               !! Spatial index in x-direction.
            integer, intent(in) :: jy
               !! Spatial index in y-direction.
            real(rk), intent(in) :: uxy(*)
               !! Block of `ns` dependent variables at spatial point `(jx,jy)`.
            real(rk), intent(out) :: rxy(*)
              !! Block `(jx,jy)` of \(R\).
         end subroutine rblock
      end interface
      real(rk), intent(inout) :: r1(*)
        !! Real work space available to this subroutine.
      real(rk), intent(in) :: rewt(*)
         !! Reciprocal error weights.
      real(rk), intent(in) :: cj
         !! Scalar used in forming the system Jacobian.
      real(rk), intent(out) :: bd(*)
         !! LU factors of the diagonal blocks.
      integer, intent(out) :: ipbd(*)
         !! Pivots for the LU factorizations.
      integer, intent(inout) :: ierr
        !! Error flag.
        !! `0`: no error occurred;
        !! `k > 0`: a zero pivot was found at the `k`-th stage in one of the LU factorizations.

      integer :: i, ibd, idiag, ig, igx, igy, iip, j, j0, j00, js, jx, jy
      real(rk) :: del, dfac, fac, uj

      ! Make MP calls to RBLOCK to approximate each diagonal block of dR/du.
      dfac = 1e-2_rk
      ibd = 0
      do igy = 1, ngy
         jy = jyr(igy)
         j00 = (jy - 1)*mxmp
         do igx = 1, ngx
            jx = jxr(igx)
            j0 = j00 + (jx - 1)*mp
            ! If R0 has not been set previously as an array of length NEQ, it can
            ! be set here, as an array of length MP, with the call
            !    CALL RBLOCK (T, JX, JY, U(J0+1), R0)
            ! In this case, change R0(J0+I) below to R0(I).
            do js = 1, mp
               j = j0 + js
               uj = u(j)
               del = max(srur*abs(uj), dfac/rewt(j))
               u(j) = u(j) + del
               fac = -one/del
               call rblock(t, jx, jy, u(j0 + 1), r1)
               do i = 1, mp
                  bd(ibd + i) = (r1(i) - r0(j0 + i))*fac
               end do
               u(j) = uj
               ibd = ibd + mp
            end do
         end do
      end do

      ! Add matrix CJ*Id, and do LU decomposition on blocks.
      ibd = 1
      iip = 1
      do ig = 1, ngrp
         idiag = ibd
         do i = 1, mp
            if (i <= mpd) bd(idiag) = bd(idiag) + cj
            idiag = idiag + (mp + 1)
         end do
         call dgefa(bd(ibd), mp, mp, ipbd(iip), ierr)
         if (ierr /= 0) exit
         ibd = ibd + mpsq
         iip = iip + mp
      end do

   end subroutine jac_rbgpre