daskr_banpre Module

Preconditioner Routines for Banded Problems

The following pair of subroutines – banja and banps – provides a general-purpose banded preconditioner matrix for use with the DASKR solver, with the Krylov linear system method. When using DASKR to solve a problem $G(t,y,y') = 0$, whose iteration matrix (Jacobian) $J = dG/dy + c_j dG/dy'$, where $c$ is a scalar, is either banded or approximately equal to a banded matrix, these routines 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 (info(12) = 1) in DASKR.

Other than the user-supplied residual routine res defining $G(t,y,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.

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 banja and banps. If the user program also uses ipar for communication with res, that data should be located beyond the first 2 positions.

• Import this module. Set info(15) = 1 to indicate that a jac routine exists. Then in the call to DASKR, pass the procedure names banja and banps as the arguments jac and psol, respectively.

• The DASKR work arrays rwork and iwork must include segments wp and iwp for use by banja and banps. The lengths of these arrays depend on the problem size and half-bandwidths, as follows:

   lwp  = length of rwork segment wp  = (2*ml + mu + 1)*neq + 2*((neq/(ml + mu + 1)) + 1)
   liwp = length of iwork segment iwp = neq

Note the integer divide in lwp. Load these lengths in iwork as iwork(27) = lwp and iwork(28) = liwp, and include these values in the declared size of rwork and iwork.

The banja and banps routines generate and solve the banded preconditioner matrix $P$ within the preconditioned Krylov algorithm used by DASKR when info(12) = 1. $P$ is generated and LU-factored periodically during the integration, and the factors are used to solve systems $Px = b$ as needed.