daskr_banpre Module

Preconditioner Routines for Banded Problems

This module provides a general-purpose banded preconditioner for use with the daskr solver, with the Krylov linear system method.

When using daskr to solve a problem , whose Jacobian , where is a scalar, is either banded or approximately equal to a banded matrix, the routines jac_banpre and psol_banpre can be used to generate a banded approximation to as the preconditioner and to solve the resulting banded linear system, in conjunction with the Krylov method option.

Other than the user-supplied residual routine res defining , the only other inputs required by these routines are the half-bandwidth parameters and of the approximate banded Jacobian. If the system size is , the half-bandwidths are defined as integers between 0 and such that only elements with indices satisfying are to be retained in the preconditioner. For example, if , a diagonal matrix will be generated as the preconditioner. The banded preconditioner is obtained by difference quotient approximations. If the true problem Jacobian is not banded but is approximately equal to a matrix that is banded, the procedure used here will have the effect of lumping the elements outside of the band onto the elements within the band.

Usage

To use these routines in conjunction with daskr, the user's calling program should include the following, in addition to setting the other daskr input parameters:

  • Dimension the array ipar to have length at least 2, and load the half-bandwidths into ipar as ipar(1) = ml and ipar(2) = mu. Array ipar is used to communicate these parameters to jac_banpre and psol_banpre. If the user program also uses ipar for communication with res, that data should be located beyond the first 2 positions.

  • Set info(12) = 1 to select the Krylov iterative method and info(15) = 1 to indicate that a jac routine exists. Then in the call to daskr, pass the procedure names jac_banpre and psol_banpre as the arguments jac and psol, respectively.

  • The daskr work arrays rwork and iwork must include segments rwp and iwp for use by jac_banpre and psol_banpre. The lengths of these arrays depend on the problem size and half-bandwidths, as as shown in the table below. Note the integer divide in lrwp. Load these lengths in iwork as iwork(27) = lrwp and iwork(28) = liwp, and include these values in the declared size of rwork and iwork, respectively.

Variable Length
lrwp (2*ml + mu + 1)*neq + 2*((neq/(ml + mu + 1)) + 1)
liwp neq

Example

The program example_heat demonstrates the use of this preconditioner.


Uses

  • module~~daskr_banpre~~UsesGraph module~daskr_banpre daskr_banpre module~daskr daskr module~daskr_banpre->module~daskr module~daskr_kinds daskr_kinds module~daskr_banpre->module~daskr_kinds module~daskr->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

Used by

  • module~~daskr_banpre~~UsedByGraph module~daskr_banpre daskr_banpre program~example_heat example_heat program~example_heat->module~daskr_banpre

Subroutines

public 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 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.

public pure subroutine psol_banpre(neq, t, y, ydot, savr, wk, cj, wght, rwp, iwp, b, epslin, ierr, rpar, ipar)

This routine solves the linear system for the banded preconditioner , given a vector , using the LU decomposition produced by jac_banpre. The solution is carried out by dgbsl from LINPACK.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: neq

Problem size.

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

Current independent variable (not used).

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

Current dependent variables (not used).

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

Current derivatives of dependent variables (not used).

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

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

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

Real work space available to this subroutine (not used).

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

Scalar used in forming the system Jacobian (not used).

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

Error weights for computing norms (not used).

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

Real work array containing the LU decomposition of .

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

Integer array containing matrix pivot information.

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

Right-hand side vector on input; solution on output.

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

Tolerance for linear system (not used).

integer, intent(inout) :: ierr

Error flag.

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

Real array used for communication between the calling program and user routines (not used).

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

Integer array used for communication between the calling program and user routines (not used).