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.
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 |
The program example_heat demonstrates the use of this preconditioner.
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. |
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.
Type | Intent | Optional | 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 |
||
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). |