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
.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
procedure(res_t) | :: | res |
User-defined residuals routine. |
|||
integer, | intent(out) | :: | ires |
Error flag set by |
||
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 |
||
real(kind=rk), | intent(inout) | :: | savr(*) |
Current residual evaluated at |
||
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: |
||
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. |
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