jac_rbgpre Subroutine

public subroutine jac_rbgpre(t, u, r0, rblock, r1, rewt, cj, bd, ipbd, ierr)

Uses

  • proc~~jac_rbgpre~~UsesGraph proc~jac_rbgpre jac_rbgpre module~linpack linpack proc~jac_rbgpre->module~linpack module~daskr_kinds daskr_kinds module~linpack->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

This routine generates and preprocesses a block-diagonal preconditioner matrix, based on the part of the Jacobian corresponding to the reaction terms of the problem, using block-grouping. It generates a matrix of the form . 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 .

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: t

Current independent variable.

real(kind=rk), intent(inout) :: u(*)

Current dependent variables.

real(kind=rk), intent(in) :: r0(*)

Current value of the vector .

private subroutine rblock(t, jx, jy, uxy, rxy)

User routine that computes a single block of .

Arguments
Type IntentOptional Attributes Name
real(kind=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(kind=rk), intent(in) :: uxy(*)

Block of ns dependent variables at spatial point (jx,jy).

real(kind=rk), intent(out) :: rxy(*)

Block (jx,jy) of .

real(kind=rk), intent(inout) :: r1(*)

Real work space available to this subroutine.

real(kind=rk), intent(in) :: rewt(*)

Reciprocal error weights.

real(kind=rk), intent(in) :: cj

Scalar used in forming the system Jacobian.

real(kind=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.


Calls

proc~~jac_rbgpre~~CallsGraph proc~jac_rbgpre jac_rbgpre proc~dgefa dgefa proc~jac_rbgpre->proc~dgefa interface~idamax idamax proc~dgefa->interface~idamax

Called by

proc~~jac_rbgpre~~CalledByGraph proc~jac_rbgpre jac_rbgpre proc~jac jac proc~jac->proc~jac_rbgpre

Source Code

   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