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 $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 either banded or approximately equal to a banded matrix, the routines pjac_banpre and psol_banpre can be used to generate a banded approximation to $J$ 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 $G(t,y,\dot{y})$, the only other inputs required by these routines are the half-bandwidth parameters $\mathrm{ml}$ and $\mathrm{mu}$ of the approximate banded Jacobian. If the system size is $\mathrm{neq}$, the half-bandwidths are defined as integers between 0 and $\mathrm{neq} -1$ such that only elements with indices $(i,j)$ satisfying $-\mathrm{ml} \le j, i \le \mathrm{mu}$ are to be retained in the preconditioner. For example, if $\mathrm{ml} = \mathrm{mu} = 0$, 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 following steps are needed:

• 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 pjac_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 pjac routine exists. Then in the call to daskr, pass the procedure names pjac_banpre and psol_banpre as the arguments pjac and psol, respectively.

• The daskr work arrays rwork and iwork must include segments rwp and iwp for use by pjac_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.

Example

The program example_heat demonstrates the use of this preconditioner.