This module provides a general-purpose sparse incomplete LU (ILU) 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 a general sparse matrix, the routines jac_ilupre and psol_ilupre can be used to generate an approximation to as the preconditioner and to solve the resulting sparse linear system, in conjunction with the Krylov method option.
Internally, the incomplete LU factorization is achieved via one of two routines — ilut
or ilutp — from the SPARSKIT library. The routine ilut performs an ILU factorization
of a sparse matrix using a dual thresholding technique based on a drop tolerance (tolilut
)
and a level of fill-in parameter (lfililut
). The parameter lfililut
controls the amount
of fill-in allowed in the factorization (limited to a maximum of 2*lfililut*neq
, but normally
much less). Increasing lfililut
will generally make the ILU factorization more accurate.
The parameter tolilut
also controls the accuracy of the ILU factorization via a drop tolerance
based on element size. Decreasing tolilut
will increase the amount of fill-in and make for
a more accurate factorization. The routine ilutp is a variant of ilut that in addition
performs pivoting based on a tolerance ratio permtol
.
An important aspect of using incomplete factorization techniques is that of reordering the
rows and columns in the Jacobian matrix before performing the ILU. In this package,
this is accomplished via the parameter ireorder
, which when equal to 1 performs a reverse
Cuthill-McKee (RCM) reordering before performing the ILU factorization. Based on the limited
amount of testing done so far, RCM seems the best overall choice. It is possible to include
a different reordering technique if desired.
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:
ipar
to have length at least 30, and load the following parameters
into ipar
as:Index | Name | Description |
---|---|---|
1 | ml |
The lower bandwidth used in calculating . |
2 | mu |
The upper bandwidth used in calculating . |
3 | lenpfac |
The average number of nonzeros in a row of . The maximum of nonzeros allowed in is nnzmx = lenpfac*neq . lenpfac >= 2 . |
4 | lenplufac |
The average amount of fill-in per row in the factored . The maximum number of nonzeros allowed in the factored is lenplumx = nnzmx + lenplufac*neq . lenplufac >=2 . |
5 | ipremeth |
Preconditioner type flag. =1 means ilut and =2 means ilutp. |
6 | lfililut |
Fill-in parameter for ilut and ilutp. The largest lfililut elements per row of the L and U factors are kept. Each row of L and U will have a maximum of lfililut elements in addition to their original number of nonzero elements. |
7 | ireorder |
Reordering flag. =0 means no reordering of rows/columns before incomplete factorization. =1 means reverse Cuthill-McKee (RCM) reordering is performed. |
8 | isrnorm |
Row norm flag. =1 means compute and use row norms as scalings in the preconditioner system . =0 means no row norm scaling. |
9 | normtype |
Type of row norm scaling for isrnorm . =0 means max-norm. =1 means 1-norm. =2 means 2-norm. |
10 | jacout |
Output Jacobian flag. =1 means write and initial residual to a file (logical unit in ipar(29) ). Integration halts with ires = -2 . =0 means no output. Storage format is Boeing-Harwell. |
11 | jscalcol |
Flag for scaling columns by the inverses of elements in the ewt array. =0 means no scaling. =1 means perform scaling. |
21:28 | — | Used to hold pointer information. |
29 | jacout_unit |
Logical unit number for matrix output file. Used only if jacout = 1 . |
30 | rescalls |
On return from daskr, holds the number of calls to the res routine used in preconditioner evaluations. |
rpar
to have length at least 2, and load the parameters shown in the
table below into rpar
.Index | Name | Description |
---|---|---|
1 | tolilut |
Drop tolerance for use by ilut and ilutp. tolilut >= 0 . Larger values cause less fill-in. Good values range from 0.001 to 0.01. |
2 | permtol |
Tolerance ratio used in determining column pivoting by ILUTP. permtol >= 0 . Good values are from 0.1 to 0.01. Two columns are permuted only if a(i,j)*permtol > a(i,i) . |
The two parameters tolilut
and lfililut
give the user a great deal of flexibility. One
can use tolilut = 0
to get a strategy based on keeping the largest elements in each row
of and . Taking tolilut /= 0
but lfililut = neq
will give the usual threshold
strategy (however, fill-in is then unpredictable).
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_ilupre
and psol_ilupre
as the arguments jac
and psol
, respectively.
The daskr work arrays rwork
and iwork
must include segments rwp
and iwp
for use
by jac_ilupre and psol_ilupre. The lengths of these depend on the problem size,
half-bandwidths, and other parameters as shown in the table below. 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*lenpfac*neq + lenplufac*neq + isrnorm*neq + neq |
liwp |
3*neq + 1 + 3*lenpfac*neq + 2*lenplufac*neq + 2*ireorder*neq + 2*(ipremeth-1)*neq |
The program example_heatilu demonstrates the use of this preconditioner.
Setup routine for the incomplete LU preconditioner. This routine checks the user input and calculates the minimum length needed for the preconditioner workspace arrays.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | neq |
Problem size. |
||
integer, | intent(in) | :: | lrwp |
Current length of |
||
integer, | intent(in) | :: | liwp |
Current length of |
||
real(kind=rk), | intent(in) | :: | rpar(*) |
User real workspace. |
||
integer, | intent(in) | :: | ipar(*) |
User integer workspace. |
||
integer, | intent(out) | :: | ierr |
Error flag (0 means success, else failure):
|
||
integer, | intent(out) | :: | lrwp_min |
Minimum |
||
integer, | intent(out) | :: | liwp_min |
Minimum |
This subroutine uses finite-differences to calculate the Jacobian matrix in sparse format, and then performs an incomplete LU decomposition using either ilut or ilutp from SPARSKIT.
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(neq) |
Current dependent variables. |
||
real(kind=rk), | intent(inout) | :: | ydot(neq) |
Current derivatives of dependent variables. |
||
real(kind=rk), | intent(in) | :: | rewt(neq) |
Reciprocal error weights for scaling |
||
real(kind=rk), | intent(inout) | :: | savr(neq) |
Current residual evaluated at |
||
real(kind=rk), | intent(inout) | :: | wk(neq) |
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(*) |
Matrix elements of ILU. |
||
integer, | intent(inout) | :: | iwp(*) |
Array indices for elements of ILU. |
||
integer, | intent(inout) | :: | ierr |
Error flag (0 means success, else failure). |
||
real(kind=rk), | intent(inout) | :: | rpar(*) |
Real array used for communication between the calling program and user routines. |
||
integer, | intent(inout) | :: | ipar(*) |
Integer array used for communication between the calling program and user routines. |
This subroutine solves the linear system for the sparse preconditioner , given a vector , using the incomplete LU decomposition produced by jac_ilupre. The solution is carried out by lusol from SPARSKIT.
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) | :: | r0(*) |
Current residual evaluated at |
||
real(kind=rk), | intent(in) | :: | wk(*) |
Real work space available to this subroutine. |
||
real(kind=rk), | intent(in) | :: | cj |
Scalar used in forming the system Jacobian. |
||
real(kind=rk), | intent(in) | :: | wght(*) |
Error weights for computing norms. |
||
real(kind=rk), | intent(inout) | :: | rwp(*) |
Matrix elements of ILU. |
||
integer, | intent(inout) | :: | iwp(*) |
Array indices for elements of ILU |
||
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 (not used). |
||
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. |