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 $G(t,y,\dot{y}) = 0$, whose Jacobian $J = \partial G/ \partial y + c_J \partial G/ \partial \dot{y}$, where $c_J$ is a scalar, is a general sparse matrix, the routines jac_ilupre and psol_ilupre can be used to generate an approximation to $J$ 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 $J$ 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:

• Dimension the array rpar to have length at least 2, and load the parameters shown in the table below into rpar.

• 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 $L$ and $U$. 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.

Example

The program example_heatilu demonstrates the use of this preconditioner.