jac_banpre Subroutine

public subroutine jac_banpre(res, ires, neq, t, y, ydot, rewt, savr, wk, h, cj, rwp, iwp, ierr, rpar, ipar)

Uses

  • proc~~jac_banpre~~UsesGraph proc~jac_banpre jac_banpre module~linpack linpack proc~jac_banpre->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 a banded preconditioner matrix that approximates the Jacobian matrix . The banded matrix has half-bandwidths and , and is computed by making 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.

Arguments

Type IntentOptional Attributes Name
procedure(res_t) :: res

User-defined residuals routine.

integer, intent(out) :: ires

Error flag set by res.

integer, intent(in) :: neq

Problem size.

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

Current independent variable.

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

Current dependent variables.

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

Current derivatives of dependent variables.

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

Reciprocal error weights for scaling y and ydot.

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

Current residual evaluated at (t, y, ydot).

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

Real work space available to this subroutine.

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

Current step size.

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

Scalar used in forming the system Jacobian.

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

Real work array for , etc. On output, it contains the LU decomposition of the banded approximation .

integer, intent(inout) :: iwp(*)

Integer work space for matrix pivot information.

integer, intent(inout) :: ierr

Error flag: ierr > 0 if is singular, and ierr = 0 otherwise.

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


Calls

proc~~jac_banpre~~CallsGraph proc~jac_banpre jac_banpre proc~dgbfa dgbfa proc~jac_banpre->proc~dgbfa interface~daxpy daxpy proc~dgbfa->interface~daxpy interface~dscal dscal proc~dgbfa->interface~dscal interface~idamax idamax proc~dgbfa->interface~idamax

Source Code

   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