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 $\langle$ method $\rangle$ (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 ${\bf N+1}$ 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 ${\bf 3*N}$ 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.⁵ [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 $\omega$ 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;
$\neq 0$:	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 ${\bf RPARM(I), I=2,3,5,6,7}$ 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

[$\neq 2$	Case I: Fixed ${\bf SME} \leq m(B)$ (general case);
=2	Case II: Use when it is known that $\vert m(B)\vert \leq M(B)]$

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 ${\bf -CME}$ 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).⁶ Also, if A is an L-matrix, then for the Jacobi matrix, we have $\vert m(B)\vert \leq M(B)$ and SME is always ${\bf -CME}$ (Case II). ⁷ 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 ${\bf RPARM(I), I=2,3}$ 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;
$\geq 0$:	skip indexing--system already in red-black form;

For other methods,

<0:	skip indexing--system already in desired form;
$\geq 0$:	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; $\neq 0$: test done]

IPARM(11)

ITIME is the timing switch. Default: 0 [0: time method; $\neq 0$: 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 ${\bf RPARM(I), I=11,12}$,
	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 ${\bf LEVEL} \leq 0$, no printing is done. See ${\bf RPARM(I), I=11,12}$ 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: $5 \times 10^{-6}$

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. ${\bf CME} \leq M(B)$. 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 $\leq m(B)$. In Case II, SME is always set to ${\bf -CME}$ 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 ${\bf IPARM(6)=1}$ 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 $\frac{1}{4}$, 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. ⁸ 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: $100.\times{\bf SRELPR}$

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()	${\bf 4*N + NCG}$
JSI()	${\bf 2*N}$
SOR()	${\bf N}$
SSORCG()	${\bf 6*N + NCG}$
SSORSI()	${\bf 5*N}$
RSCG()	${\bf N + 3*NB + NCG}$
RSSI()	${\bf N + NB}$

The value of NCG is ${\bf 2*IPARM(1)}$ for symmetric sparse storage and ${\bf 4*IPARM(1)}$ 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
${\bf IER} =$	0,	Normal convergence was obtained.
=	$1+{\rm Mth}$,	Invalid order of the system, N.
=	$2+{\rm Mth}$,	Workspace array WKSP(*) is not large enough. IPARM(8)
		is set to the amount of required workspace, NW.
=	$3+{\rm Mth}$,	Failure to converge in IPARM(1) iterations. RPARM(1)
		is reset to the last stopping value computed.
=	$4+{\rm Mth}$,	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.⁹