Next: User-Oriented Modules
Up: ITPACK 2C: A FORTRAN
Previous: Sparse Matrix Storage
Usage
The user is expected to provide the coefficient matrix and the
right-hand side of the linear system to be solved. The data structure
for the matrix of the system is either the symmetric or nonsymmetric
sparse storage format described in Section 2. An initial
guess for the solution should be provided, if one is known; otherwise,
it can be set to all zero values. A series of approximations for the
solution are generated iteratively until the convergence criteria is
satisfied. The algorithms are performed in two work space arrays and
some control over the algorithmic procedure can be obtained from
switches in two parameter arrays.
There are seven main subroutines in ITPACK, each corresponding to an
iterative method. They are:
Subroutine |
Method |
JCG() |
Jacobi Conjugate Gradient |
JSI() |
Jacobi Semi-Iteration |
SOR() |
Successive Overrelaxation |
SSORCG() |
Symmetric SOR Conjugate Gradient |
SSORSI() |
Symmetric SOR Semi-Iteration |
RSCG() |
Reduced System Conjugate Gradient |
RSSI() |
Reduced System Semi-Iteration |
and the calling sequence is:
CALL
method
(N, IA, JA, A, RHS, U, IWKSP,
NW, WKSP, IPARM,
RPARM, IER)
where the parameters are defined in the following. Here ``input" means
that the subroutine expects the user to provide the necessary input data
and ``output" means that the routine passes back information in the
variable or array indicated. All parameters are linear arrays except
variables N, NW, and IER. Moreover, all parameters may
be altered by the subroutine call except variables N and
NW. (See Section 7 for additional details.)
- N
- is the order of the linear system. [integer; input]
- IA(*)
- is a vector of length
used in the sparse
matrix storage format. It contains the row pointers
into JA(*) and A(*). [integer array; input]
- JA(*)
- is a vector of length NZ (defined in A(*) below)
used in the sparse matrix storage format. It contains
the column numbers for the corresponding entries in
A(*). [integer array; input]
- A(*)
- is a vector of length NZ used in the sparse matrix
storage format. It contains the nonzero entries of the
coefficient matrix with positive diagonal elements.
(NZ is the number of nonzero entries in the upper
triangular part of the coefficient matrix when symmetric
storage is used and is the total number of nonzeros when
nonsymmetric storage is used.) [real array; input]
- RHS(*)
- is a vector of length N containing the right-hand side
of the linear system. [real array; input]
- U(*)
- is a vector of length N containing the initial guess
to the solution of the linear system on input and the
latest approximate solution on output. [real array;
input/output]
- IWKSP(*)
- is a vector of length
used for integer workspace.
When reindexing for red-black ordering, the first
N locations contain on output the permutation
vector for the red-black indexing, the next N
locations contain its inverse, and the last N
are used for integer workspace.5
[integer array; output]
- NW
- is a scalar. On input, NW is the available length for
WKSP(*). On output, IPARM(8) is the actual
amount used (or needed). [integer; input]
- WKSP(*)
- is a vector used for real working space whose length depends
on the iterative method being used. It must be at least
NW entries long. (See the table near the end of this
section for the required amount of workspace for each
method.) [real array]
- IPARM(*)
- is a vector of length 12 used to initialize various
integer and logical parameters. Default values may be
set by calling subroutine DFAULT() described below.
On output, IPARM(*) contains the values of the
parameters that were changed. (Further details are
given later in this section.) [integer array;
input/output]
- RPARM(*)
- is a vector of length 12 used to initialize various real
parameters on input. Default values may be set by
calling subroutine DFAULT() described below. On
output, RPARM(*) contains the final values of
the parameters that were changed. (Further details
are given later in this section.) [real array;
input/output]
- IER
- is the error flag which is set to zero for normal convergence
and to a nonzero integer when an error condition is
present. (See the table at the end of this section
for the meaning of nonzero values.) [integer; output]
The user may supply nondefault values for selected quantities in
IPARM(*) and
RPARM(*) by first executing
CALL DFAULT (IPARM, RPARM)
and then assigning the appropriate nondefault values before calling a
solution module of ITPACK.
The iterative algorithms used in ITPACK are quite complicated and some
knowledge of iterative methods is necessary to completely understand
them. The interested reader should consult the technical report [4]
and the book [6] for details. Important variables in this package
which may change adaptively are CME (estimate of M(B), the largest
eigenvalue of the Jacobi matrix), SME (estimate of m(B), the
smallest eigenvalue of the Jacobi matrix), OMEGA (overrelaxation
parameter for the SOR and SSOR methods), SPECR (estimated spectral
radius of the SSOR matrix), BETAB (estimate for the spectral radius
of the matrix LU where L and U are strictly lower and upper
triangular matrices, respectively, such that the Jacobi matrix B=L+U).
The integer array IPARM(*) and real array RPARM(*) allow
the user to control certain parameters which affect the performance
of the iterative algorithms. Furthermore, these arrays allow the updated
parameters from the automatic adaptive procedures to be communicated
back to the user. The entries in IPARM(*) and RPARM(*) are:
- IPARM(1)
- ITMAX is the maximum number of iterations
allowed. It is reset on output to the number of
iterations performed. Default: 100
- IPARM(2)
- LEVEL is used to control the level of output.
Each higher value provides additional information.
Default: 0
[<0: |
no output on unit IPARM(4); |
0: |
fatal error messages only; |
1: |
warning messages and minimum output; |
2: |
reasonable summary (progress of algorithm); |
3: |
parameter values and informative comments; |
4: |
approximate solution after each iteration
(primarily useful for debugging); |
5: |
original system] |
- IPARM(3)
- IRESET is the communication switch. Default: 0
[0: |
implies certain values of IPARM(*) and
RPARM(*) will be overwritten |
|
to communicate parameters back to the user; |
: |
only IPARM(1) and IPARM(8)
will be reset.] |
- IPARM(4)
- NOUT is the output unit number. Default: 6
- IPARM(5)
- ISYM is the sparse storage format switch.
Default: 0
[0: |
symmetric sparse storage; |
1: |
nonsymmetric sparse storage] |
- IPARM(6)
- IADAPT is the adaptive switch. It determines
whether certain parameters have been set by the user or
should be computed automatically in either a fully
or partially adaptive sense. Default: 1
[0: |
fixed iterative parameters used for SME,
CME, OMEGA, SPECR, and |
|
BETAB (nonadaptive); |
1: |
fully adaptive procedures used for all parameters; |
2: |
(SSOR methods only) SPECR determined
adaptively and CME, BETAB, |
|
and OMEGA fixed; |
3: |
(SSOR methods only) BETAB fixed and all other
parameters determined |
|
adaptively] |
(See [4,6] for details and
for CME, SME, OMEGA, SPECR,
BETAB, respectively. These parameters are set
by subroutine DFAULT() or by the user.)
- IPARM(7)
- ICASE is the adaptive procedure case switch for the
JSI and SSOR methods. There are two strategies, called
Case I and Case II, for doing the adaptive procedure.
The choice of which case to select corresponds to knowledge
of the eigenvalues of the Jacobi matrix B and their
estimates. Default: 1
[ |
Case I: Fixed
(general case); |
=2 |
Case II: Use when it is known that
|
The case switch determines how the estimates for
SME and CME are recomputed adaptively. In
Case I, SME is fixed throughout and should be less
than or equal to m(B). In Case II, SME is set
to
which may adaptively change. As far as
the adaptive procedure is concerned, Case I is the most
general case and should be specified in the absence of
specific knowledge of the relationship between the
eigenvalues and their estimates. An example when Case II
is appropriate occurs when the Jacobi matrix has
Property A, since
m(B) = -M(B).6 Also, if A is an L-matrix,
then for the Jacobi matrix, we have
and SME is always
(Case II).
7 Selecting the
correct case may increase the rate of convergence of
the iterative method. (See [6] for additional
discussion on Case I and II. Also, see
for CME, SME,
respectively.)
- IPARM(8)
- NWKSP is the amount of workspace used. It is used
for output only. If ITMAX is set to a value just
over the actual number of iterations necessary for
convergence, the amount of memory for WKSP(*) can
be reduced to just over the value returned here. This
may be done when rerunning a problem, for example.
Default: 0
- IPARM(9)
- NB is the red-black ordering switch. On output,
if reindexing is done, NB is set to the order of
the black subsystem. Default: -1
[For RS methods,
<0: |
compute red-black indexing and permute
system; |
: |
skip indexing--system already in red-black
form; |
For other methods,
<0: |
skip indexing--system already in desired
form; |
: |
compute red-black indexing and permute
system] |
A negative integer value for IPARM(9) causes the
equations to be handled in the most general way
appropriate for the solution method being used. For
methods other than RS methods this is the ``natural
order" while for RS methods it is the ``red-black order."
A non-negative value produces a red-black permutation
for all methods except for the RS methods which are
assumed to be in red-black order with the order of the
black subsystem NB given. If reindexing is performed,
IPARM(9) will contain the order of the black subsystem
on output.
- IPARM(10)
- IREMOVE is the switch for effectively removing
rows and columns when the diagonal entry is extremely
large compared to the nonzero off-diagonal entries in
that row. (See RPARM(8) for additional details.)
Default: 0
[0: not done; :
test done]
- IPARM(11)
- ITIME is the timing switch. Default: 0
[0: time method; :
not done]
- IPARM(12)
- IDGTS is the error analysis switch. It determines
if an analysis is done to determine the accuracy of
the final computed solution. Default: 0
[<0: |
skip error analysis |
0: |
compute DIGIT1 and DIGIT2 and
store in
, |
|
respectively; |
1: |
print DIGIT1 and DIGIT2; |
2: |
print final approximate solution vector; |
3: |
print final approximate residual vector; |
4: |
print both solution and residual vectors; |
|
otherwise: no printing] |
(If
,
no printing is done. See
for details on DIGIT1
and DIGIT2.)
- RPARM(1)
- ZETA is the stopping criterion or approximate
relative accuracy desired in the final computed solution.
If the method did not converge in IPARM(1)
iterations, RPARM(1) is reset to an estimate of
the relative accuracy achieved. The stopping criterion
is a test of whether ZETA is greater than the
ratio of the two norm of the pseudo-residual vector and
the two norm of the current iteration vector times a
constant involving an eigenvalue estimate. (See [4,6]
for details.) Default:
- RPARM(2)
- CME is the estimate of the largest eigenvalue of
the Jacobi matrix. It changes to a new estimate if the
adaptive procedure is used.
.
Default: 0.0
- RPARM(3)
- SME is the estimate of the smallest eigenvalue of
the Jacobi matrix for the JSI method. In Case I,
SME is fixed throughout at a value .
In Case II, SME is always set to
with CME changing in the adaptive procedure. (See
IPARM(7) for definitions of Case I and II.)
Default: 0.0
- RPARM(4)
- FF is the adaptive procedure damping factor. Its
values are in the interval (0.,1.] with 1. causing
the most frequent parameter changes when the fully
adaptive switch
is used.
Default: 0.75
- RPARM(5)
- OMEGA is the overrelaxation parameter for the SOR
and the SSOR methods. If the method is fully adaptive,
OMEGA changes. Default: 1.0
- RPARM(6)
- SPECR is the estimated spectral radius for the
SSOR matrix. If the method is adaptive, SPECR
changes. Default: 0.0
- RPARM(7)
- BETAB is the estimate for the spectral radius of the
matrix LU used in the SSOR methods. BETAB may
change depending on the adaptive switch IPARM(6).
The matrix L is the strictly lower triangular part of
the Jacobi matrix and U is the strictly upper triangular
part. When the spectral radius of LU is less than or
equal to
,
the ``SSOR condition" is satisfied
for some problems provided one uses the natural ordering.
(See [4,5,18] for additional details.) Default: 0.25.
- RPARM(8)
- TOL is the tolerance factor near machine relative
precision, SRELPR. In each row, if all nonzero
off-diagonal row entries are less than TOL times
the value of the diagonal entry, then this row and
corresponding column are essentially removed from the
system. This is done by setting the nonzero off-diagonal
elements in the row and corresponding column to zero,
replacing the diagonal element with 1, and adjusting
the elements on the right-hand side of the system so that
the new system is equivalent to to the original one.
8
If the diagonal entry is the only nonzero
element in a row and is not greater than the reciprocal
of TOL, then no elimination is done. This procedure
is useful for linear systems arising from finite element
discretizations of PDEs in which Dirichlet boundary
conditions are handled by giving the diagonal values in
the linear system extremely large values. (The installer
of this package should set the value of SRELPR.
See comments in subroutine DFAULT() for additional
details.) Default:
- RPARM(9)
- TIME1 is the total time in seconds from the beginning
of the iterative algorithm until convergence. (A machine
dependent subprogram call for returning the time in
seconds is provided by the installer of this package.)
Default: 0.0
- RPARM(10)
- TIME2 is the total time in seconds for the entire
call. Default: 0.0
- RPARM(11)
- DIGIT1 is the approximate number of digits using
the estimated relative error with the final approximate
solution. It is computed as the negative of the logarithm
base ten of the final value of the stopping test. (See
details below or [6].) Default: 0.0
- RPARM(12)
- DIGIT2 is the approximate number of digits using
the estimated relative residual with the final approximate
solution. It is computed as the negative of the logarithm
base ten of the ratio of the two norm of the residual
vector and the two norm of the right-hand side vector.
This estimate is related to the condition number of the
original linear system and, therefore, it will not be
accurate if the system is ill-conditioned. (See details
below or [6].) Default: 0.0
DIGIT1 is determined from the actual stopping test computed on
the final iteration, whereas DIGIT2 is based on the computed
residual vector using the final approximate solution after the algorithm
has converged. If these values differ greatly, then either the stopping
test has not worked successfully or the original system is ill-conditioned.
(See [6] for additional details.)
For storage of certain intermediate results, the solution modules
require a real vector WKSP(*) and a corresponding variable NW
indicating the available space. The length of the workspace array varies
with each solution module and the maximum amount needed is given in
the following table.
Solution Module |
Maximum Length of WKSP(*) |
JCG() |
|
JSI() |
|
SOR() |
|
SSORCG() |
|
SSORSI() |
|
RSCG() |
|
RSSI() |
|
The value of NCG is
for symmetric sparse storage
and
for nonsymmetric sparse storage. It should be noted
that the actual amount of workspace used may be somewhat less than these
upper limits since some of the latter are dependent on the maximum number
of iterations allowed, ITMAX, stored in IPARM(1). Clearly,
the array WKSP(*) must be dimensioned to at least the value of
NW.
Nonzero integer values of the error flag IER indicate that an error
condition was detected. These values are listed below according to
their numerical value and to the name of the routine in which the flag
was set.
Error Flag |
Meaning |
|
0, |
Normal convergence was obtained. |
= |
, |
Invalid order of the system, N. |
= |
, |
Workspace array WKSP(*) is not large
enough. IPARM(8) |
|
|
is set to the amount of required workspace,
NW. |
= |
, |
Failure to converge in IPARM(1)
iterations. RPARM(1) |
|
|
is reset to the last stopping value computed. |
= |
, |
Invalid order of the black subsystem, NB. |
= |
101, |
A diagonal element is not positive. |
= |
102, |
No diagonal element in a row. |
= |
201, |
Red-black indexing is not possible. |
= |
301, |
No entry in a row of the original matrix. |
= |
302, |
No entry in a row of the permuted matrix. |
= |
303, |
Sorting error in a row of the permuted matrix. |
= |
401, |
A diagonal element is not positive. |
= |
402, |
No diagonal element in a row. |
= |
501, |
Failure to converge in ITMAX
function evaluations. |
= |
502, |
Function does not change sign at the endpoints. |
= |
601, |
Successive iterates are not monotone increasing. |
JCG(), JSI(), SOR(), SSORCG(), SSORSI(),
RSCG(), RSSI() assign values to Mth of 10,20,30,40,50,60,70,
respectively. SBELM(), PRBNDX(), PERMAT(),
SCAL(), ZBRENT(), EQRT1S() are subroutines with
error flags in the 100's, 200's, 300's, 400's, 500's, 600's, respectively.
These routines perform the following functions: SBELM() removes
rows and columns, PRBNDX() determines the red-black indexing,
SCAL() scales the system, ZBRENT() is a modified IMSL routine
for computing a zero of a function which changes sign in a given
interval, EQRT1S() is a modified IMSL routine for computing the
largest eigenvalue of a symmetric tridiagonal matrix.9
Next: User-Oriented Modules
Up: ITPACK 2C: A FORTRAN
Previous: Sparse Matrix Storage