parmdesc.inc
c* Description of the parameter arrays iparm and rparm. Defaults
c* set by a call to JT_SetDefaults are indicated.
c*
c* iparm(1): [not yet implemented]
c* safety - toggles extra input checking
c* 0 ==> no extra checking
c* 1 ==> extra checking
c*
c* iparm(2):
c* iter - number of iterations completed
c* Note that if iter=0 and status=0, the initial guess was
c* found to be the solution, so no iterations were required.
c*
c* iparm(3):
c* out - output mode:
c* 0 ==> no output
c* 1 ==> write out error messages only
c* 2 ==> write out error and warning messages [DEFAULT]
c* 3 ==> same as iout=2, plus write a one-line summary
c* consisting of the iteration number, the norms of
c* the calculated residual and error estimate, and the
c* norms of the true residual and error estimate (if
c* calculated) each time convergence is checked
c* 4 ==> same as iout=3, plus write out the coefficient,
c* preconditioner (if there is one), source, initial
c* guess and converged solution
c* 5 ==> same as iout=4, plus write out the current iterate
c* and true residual (if computed) at each iteration
c* 6 ==> same as iout=5, plus write out vectors at each stage
c* of the iteration, details about the preconditioner
c* if available, etc. (WARNING: this generates reams of
c* output)
c*
c* iparm(4):
c* luout - logical unit to which output should be written
c* [DEFAULT is 6, which is typically stdout]
c*
c* iparm(5):
c* luerr - logical unit to which error messages should be
c* written [DEFAULT should be stderr, so is set to 0 except
c* on HP, which uses unit 7 for stderr]
c*
c* iparm(6):
c* itmax - maximum number of iterations to allow [NO DEFAULT]
c*
c* NOTE: itmax must be less than or equal to the maximum dimension
c* of rnorm, err, rnormt, and errt, and must be set.
c*
c* iparm(7):
c* stop - determines which stopping test to use
c* (if negative, then the true residual will be used)
c* 0 ==> || x - xold || / || x || [DEFAULT]
c* 1 ==> || r || / ( || A ||*|| x || + || b || )
c* 2 ==> || r || / || b ||
c* 3 ==> || r || / || x ||
c* 4 ==> || r || / || r0 ||
c* 5 ==> || r ||
c*
c* where A, x, and b are the coefficient matrix, current estimate
c* of the solution vector, and source vector of the linear system
c* being solved, Ax = b. The vector r is the residual, r = b - Ax,
c* and r0 is the initial residual (the residual computed using the
c* initial guess for x).
c*
c* Recommendations: CG -> stop = 2
c* GMRES -> stop = 2
c* BCGS -> stop = -2 (see note)
c* TFQMR -> stop = -2 (see note)
c* Jacobi -> stop = 0
c* SOR -> stop = 0
c* SSOR -> stop = 0
c*
c* NOTES:
c* - The norm to use is set by iparm(8).
c*
c* - Stopping test #0 is fine for the stationary methods (Jacobi,
c* SOR, SSOR) since 1) the residual is not directly available,
c* and 2) these methods converge nearly monotonically anyway.
c* However, this test is completely inappropriate and dangerous
c* for the nonstationary methods. It is the default only
c* because it is the only test that "works" without extra work
c* for all methods.
c*
c* - Stopping test #2 is probably the best for most situations.
c* However, it can lead to difficulties when
c* || A ||*|| x || >> || b ||
c* In that case, use stopping test #1.
c*
c* - Stopping tests #1 and #3 increase the cost of GMRES
c* significantly, since a new iterate is normally not computed
c* at each iteration.
c*
c* - Note that if stopping test #1 is used, the norm of the
c* coefficient is needed. This value is calculated unless the
c* storage format is unknown (storage = 0), in which case
c* the value must be passed in as rparm(4).
c*
c* - Stopping test #4 is not recommended, as it depends too much
c* on the initial guess. Also note that when the initial guess
c* of x = 0 is used, this test is equivalent to test #2.
c*
c* - Stopping test #5 should be used carefully (if ever).
c*
c* - Stopping test #2, based on the true residual, is recommended
c* for BCGS and TFQMR simply to protect against false
c* convergence. However, it adds significant overhead to each
c* iteration.
c*
c* A better (more efficient) way to protect against false
c* convergence is to use a stopping test based on the
c* recursively-generated residual (stop > 0), and then compute
c* the true residual and corresponding error estimate after
c* "convergence". If the error estimate based on the true
c* residual does not meet the desired convergence criterion, the
c* the final iterate can be used as the initial guess for
c* additional iterations with the same (or a different) method.
c*
c* - For the stationary iterative methods (Jacobi, SOR, SSOR), the
c* true residual is always used when stop /= 0, since an
c* estimated residual isn't built up during the iteration. So,
c* e.g., stop = 1 is the same as stop = -1. This adds quite a
c* bit of expense to each iteration, but is useful if you want
c* to plot the residual history.
c*
c* iparm(8):
c* norm - determines which norm to use
c* 0 ==> infinity norm
c* 1 ==> 1-norm
c* 2 ==> 2-norm for vectors, Frobenius norm for matrices
c* [DEFAULT]
c*
c* NOTES:
c* - The 2-norm is the only norm currently implemented for all
c* matrix storage formats.
c*
c* - For a vector x of length n:
c*
c* infinity norm ==> max | x |
c* 1<=i<=n i
c*
c* ---
c* \
c* 1-norm ==> / | x |
c* --- i
c* n
c*
c* ---
c* \ 2 T
c* 2-norm ==> sqrt / x = sqrt (x x)
c* --- i
c* n
c*
c* - For an n-by-m matrix A:
c* m
c* ---
c* \
c* infinity norm ==> max / | a |
c* 1<=i<=n --- ij
c* j=1
c* n
c* ---
c* \
c* 1-norm ==> max / | a |
c* 1<=j<=m --- ij
c* i=1
c* n m
c* --- ---
c* \ \ 2
c* Frobenius norm ==> sqrt / / a
c* --- --- ij
c* i=1 j=1
c*
c* - Note that the matrix 1-norm is expensive to calculate for
c* matrices stored in some formats (e.g. ELL).
c*
c* iparm(9):
c* resid - determines how often to compute and/or use true
c* residual (b-Ax) and how often to compute error estimate
c* based on true residual (errt)
c* 0 ==> never [DEFAULT]
c* 1 ==> compute b-Ax and errt each time a new iterate is
c* obtained and use as appropriate for the method
c* 2 ==> replace r with b-Ax every MAX(10,SQRT(n)) iterations
c* 3 ==> compute b-Ax every iteration, but only replace r
c* with b-Ax every MAX(10,SQRT(n)) iterations
c*
c* NOTES:
c* - Resid = 1 and resid = 3 are mainly for investigating
c* the behavior of solvers, since time is wasted calculating
c* but not using b-Ax.
c*
c* - Use of resid > 1 can be very beneficial to CG, and is
c* recommended. For production use and when CPU time is a
c* consideration, use resid = 2.
c*
c* - Values of resid > 1 are ignored by GMRES, since b-Ax
c* is recalculated automatically each nold iterations.
c*
c* - Values of resid > 1 are meaningless for stationary
c* methods (Jacobi, SOR, SSOR), since they do not generate
c* a residual during the course of the algorithm.
c*
c* - Values of resid > 1 are ignored by BCGS and TFQMR.
c* It doesn't seem to benefit these methods, and seems to harm
c* both in some situations.
c*
c* iparm(10):
c* method - determines which method to use if multiple methods
c* are implemented in one routine [DEFAULT = 0]
c* For JT_CG: unused
c* For JT_BCGS: 0 ==> BCGS
c* 1 ==> TFQMR
c* For JT_GMRES: unused
c* For JT_Stationary: -1 ==> Jacobi
c* 0 ==> SOR
c* 1 ==> SSOR
c*
c* NOTE: There's really not much point in using BCGS, except for
c* research purposes, since TFQMR is more robust and not much
c* more expensive computationally.
c*
c* iparm(11):
c* prein - indicates whether or not a preconditioner will be
c* passed in
c* 0 ==> no [DEFAULT]
c* 1 ==> yes (can be useful in, for example, time-dependent
c* problems when many systems with similar coefficient
c* matrices are solved)
c*
c* iparm(12):
c* pre - controls preconditioning
c* -1 ==> reverse communication (in this case control returns
c* to the calling routine whenever preconditioning is
c* required)
c* 0 ==> none [DEFAULT]
c* 1 ==> m-step Jacobi iteration
c* (also uses istep and omega)
c* 2 ==> m-step SSOR iteration
c* (also uses istep and omega)
c* 3 ==> Incomplete Cholesky
c* (also uses omega)
c* 4 ==> Incomplete LU
c* (also uses omega)
c*
c* NOTE: The incomplete factorization preconditioners are currently
c* only implemented for the full and ELL storage formats
c* (storage values of 1 and 8, respectively).
c*
c* Recommendations:
c* - M-step Jacobi (ipre=3) vectorizes well, but is typically not
c* very effective. When istep=1 it is equivalent to simply
c* diagonally scaling the original system prior to solution.
c*
c* - M-step SSOR is often very effective, but is difficult to
c* vectorize (the implementation in JTpack77 is not
c* vectorized).
c*
c* - For the incomplete factorization preconditioners, omega
c* controls the degree to which fill-in elements are added back
c* into the diagonal:
c* omega = 0 ==> standard incomplete factorization
c* omega = 1 ==> modified incomplete factorization
c* 0<omega<1 ==> relaxed incomplete factorization
c*
c* - The incomplete factorization preconditioners can be quite
c* effective in reducing the number of iterations, but are
c* often so expensive to compute as to negate any benefit in
c* using them. In addition both the generation and application
c* of incomplete factorization preconditioners is difficult to
c* vectorize, and the implementations in JTpack77 are not
c* vectorized. One way to make them more useful when solving a
c* number of linear systems with similar coefficients (e.g. a
c* time-dependent problem) is to amortize the cost of computing
c* them over several systems. That is, call the solver routine
c* initially with iprein=0, and have the solver generate and
c* use one of the incomplete factorization preconditioners.
c* Then on subsequent calls to the solver, set iprein=1, so
c* that the already-computed preconditioner will be reused.
c*
c* - While incomplete Cholesky (IC) factorization preconditioners
c* are appropriate only for symmetric coefficients, incomplete
c* LU (ILU) factorization preconditioners can be used for
c* either symmetric or nonsymmetric coefficients. Typically,
c* for a symmetric coefficient, one of the IC factorization
c* should be preferrable to one of the ILU factorization, since
c* the number of operations required to generate an IC
c* factorization would be about 1/2 that required for an ILU
c* factorization. However, if the coefficient is stored in ELL
c* format, the nature of the data structure causes generation
c* of an IC factorization to be highly inefficient. Hence, for
c* coefficients stored in ELL format, ILU factorization
c* preconditioners are preferrable to IC factorization
c* preconditioners, even if the coefficient is symmetric.
c*
c* iparm(13):
c* step - number of Jacobi or SSOR steps to take when ipre = 1
c* or 2 [DEFAULT = 1]
c*
c* iparm(14):
c* nold - number of old vectors to use in truncated methods like
c* GMRES (not required for non-truncated methods)
c* [DEFAULT = 20]
c*
c* NOTE: Typically, the number of iterations (and the likelihood of
c* failure to converge) decreases as nold increases.
c* However, the cost of the algorithm, both in operations and
c* memory required, increases as nold increases.
c*
c* rparm(1):
c* eps - convergence criterion [DEFAULT = 1e-7]
c*
c* rparm(2):
c* epspre - convergence criterion for solution of preconditioning
c* equations [DEFAULT = 1e-6]
c*
c* rparm(3):
c* omega - relaxation parameter [DEFAULT = 1.0]
c*
c* rparm(4):
c* anorm - appropriate norm of coefficient matrix
c*
c* NOTES:
c* - This parameter is only used when stop = +/- 1
c*
c* - If stop = +/- 1, the value is calculated unless the
c* storage format is unknown (storage = 0), in which case
c* the value must be passed in.
c*
c* rparm(5):
c* tiny - smallest floating-point value to worry about
c* [DEFAULT = 0.0]
c*
c* NOTES:
c* - If the initial residual norm is smaller than this value,
c* convergence is assumed.