JTPACK77 allows the host code to perform preconditioning via reverse communication. That is, the routine saves its state, sets status to some value, and returns. The calling routine thus needs to check the value of status and depending on the value, apply the preconditioner, continue, print an error message, or whatever.
To use this option, first write a routine to apply the preconditioner (with whatever arguments you wish, common blocks, or whatever). Here's an example:
#include "iadefines.h"
c**********************************************************************
c*
c* Applies preconditioner. That is, it essentially solves the
c* system:
c* Ap x = b
c*
c* <PARAMETER LIST>
c*
c* Input:
c* a - coefficient matrix
c* ia - integer vector containing info about how "a" is stored
c* NOTE: ia is described in src/include/iadesc.inc
c* ja - integer array that maps columns in a to true columns
c* ap - preconditioning matrix
c* iap - integer vector containing info about how "ap" is stored
c* NOTE: elements of iap are analogous to those in ia
c* jap - column map for preconditioning matrix
c* b - source vector
c*
c* In/Out:
c* x - vector of unknowns
c*
c* Output:
c* status - return status
c*
c**********************************************************************
subroutine Precondition (a, ia, ja, ap, iap, jap, b, x, status)
implicit none
c
c ... Input:
integer ia(*), iap(*)
integer ja(ia(_JT_leading_dimension_),*)
integer jap(ia(_JT_leading_dimension_),*)
real a(ia(_JT_leading_dimension_),*)
real ap(iap(_JT_leading_dimension_),*), b(*)
real b(*)
c
c ... In/Out:
real x(*)
c
c ... Output:
integer status
[code to perform preconditioning]
return
end
Note that in this example I have used cpp to include the file iadefines.h, which allows me to use the cpp variable _JT_leading_dimension_ for that element of the ia and iap arrays, rather than the explicit value, currently 3.
Also, I'm using those elements as the leading dimensions for arrays a, ja, ap, and jap. This is nonstandard Fortran 77 (but legal Fortran 90), and although many F77 compilers allow it (e.g.Sun Fortran 1.4 and 2.0.1, Cray Fortran), some do not (e.g.Sun Fortran 3.0.1 and later). If your Fortran compiler is strict, you'll have to dimension these arrays to (*), and pass them, along with ia(_JT_ leading_dimension_) , to a subroutine. Inside that subroutine they can be dimensioned and accessed with the desired shape.
Next, add code to the routine that calls the JTPACK77 solver routine(s) to check the return status and apply preconditioning as appropriate. Here's a chunk of code that might be used to use reverse communication with solver routine JT_ BCGS:
20 call JT_BCGS (b, a, ia, ja, iparm, rparm, x,
& cpu, rnormt, errt, rnorm, err, ap, iap, jap, status)
if (status .gt. 0) then
call Precondition (a, ia, ja, ap, iap, jap, b, x, status)
[make sure to put the result of preconditioning in x]
goto 20
elseif (status .lt. 0) then
[code to check return status for failure, as in the example using JT_ GMRES above]
endif
Note that the source vector for the preconditioning system is returned in vector b and the unknown should be returned to JT_ BCGS in vector x.
Finally, to use reverse communication for preconditioning, set iparm(_JT_ pre_)) = _JT_ pre_ reverse_ prior to calling the solver routine. Using current values of the cpp variables, this is equivalent to setting iparm(12) = -1.