daskr_ilupre Module

Preconditioner Routines for Sparse Problems

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.

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 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.
  • Dimension the array 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

Example

The program example_heatilu demonstrates the use of this preconditioner.


Uses

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

Used by

  • module~~daskr_ilupre~~UsedByGraph module~daskr_ilupre daskr_ilupre program~example_heatilu example_heatilu program~example_heatilu->module~daskr_ilupre

Subroutines

public pure subroutine setup_ilupre(neq, lrwp, liwp, rpar, ipar, ierr, lrwp_min, liwp_min)

Setup routine for the incomplete LU preconditioner. This routine checks the user input and calculates the minimum length needed for the preconditioner workspace arrays.

Arguments

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

Problem size.

integer, intent(in) :: lrwp

Current length of rwp.

integer, intent(in) :: liwp

Current length of iwp.

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): 1 <= ierr <= 11 means there's an illegal value for ipar(ierr); ierr = 12 means ipar(29) is illegal; ierr = 21 means rpar(1) is illegal; ierr = 22 means rpar(2) is illegal; ierr = 30 means more wp length is needed; ierr = 31 means more iwp length is needed.

integer, intent(out) :: lrwp_min

Minimum rwp length needed.

integer, intent(out) :: liwp_min

Minimum iwp length needed.

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

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.

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(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 y and ydot.

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

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

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.

public subroutine psol_ilupre(neq, t, y, ydot, r0, wk, cj, wght, rwp, iwp, b, epslin, ierr, rpar, ipar)

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.

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) :: r0(*)

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

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.