13 - C H A P T E R -

C H A P T E R 13

Sparse Linear System Solvers

Sun S3L provides support for solving sparse linear systems of the type A*x = B. The discussion of these routines is organized into the following sections:

Solving Sparse Linear Systems by the Direct Method

Solving Sparse Linear Systems by an Iterative Method

Deallocating a Sparse Linear System Solver

Solving Sparse Linear Systems by the Direct Method

S3L_sparse_solve solves a linear system of equations A*x = y, where A is a sparse Sun S3L array and A and y are both single- or double-precision real parallel arrays.

When calling S3L_sparse_solve to solve a new (unfactored) sparse linear system, specify S3L_FULL_FACTOR_SOLVE as the first element of the option argument vector. This will cause S3L_sparse_solve to reduce fill by reordering the array and to perform symbolic and numeric factoring before solving the system. It will also return a setup value that identifies the internal setup created by the factoring process.

If the same linear system is to be solved again, but with a different right-hand side, specify S3L_SOLVE_ONLY as the first element of the option argument. Also specify the setup value returned by the S3L_sparse_solve call that factored the sparse array. The new solution will make use of the internal setup created by the earlier S3L_sparse_solve call.

If a previously factored sparse array contains new values, but the sparsity pattern has not changed, it can be solved without specifying S3L_FULL_FACTOR_SOLVE. Instead, specify S3L_SAME_SPARSITY_SOLVE and the previously returned setup value. This causes S3L_sparse_solve to perform numeric factorization on the sparse array and then solve the linear system.

When the internal setup for a linear system is no longer needed, the resources associated with it can be freed by calling S3L_sparse_solve_free and specifying the applicable setup value.

Note - The S+ message-passing direct sparse solver was developed by Kai Shen and Tao Yang of the University of California at Santa Barbara. S+ can be used for general (asymmetric) sparse matrices.

The Sun Performance Library direct solver solves a sparse linear system on a single process.

S3L_sparse_solve has the following argument syntax:

S3L_sparse_solve(A, y, options, roptions, setup, ier)

A is a Sun S3L array handle that describes a parallel, single- or double-precision, real sparse array.

y is a Sun S3L array handle that describes a parallel real array of rank 1 (a vector) or rank 2 (a matrix), which contains the right-hand side of the linear system A*x = y. On exit, y is overwritten with the solution of the system.

options is an array of rank 1 whose elements control S3L_sparse_solve behavior in the following ways:

options[0] = S3L_FULL_FACTOR_SOLVE

Perform fill-reducing reordering, symbolic factorization, numeric factorization, and solve the linear system.

options[0] = S3L_SAME_SPARSITY_SOLVE

Do only numeric factorization and solve the linear system. Use this option when the sparsity pattern of a previously factored array stays the same but has a new set of values.

options[0] = S3L_SOLVE_ONLY

Only solve the linear system, using a previously computed factorization.

options[1] = S3L_SPLUS_SOLVER

Use S+ sparse solver.

options[1] = S3L_PERFLIB_SOLVER

Use the Sun Performance Library 6.0 sparse solver.

If options[1] = S3L_PERFLIB_SOLVER, specify the following options as well:

options[2] = S3L_NON_SYMMETRIC

The sparse array has asymmetric structure and asymmetric values.

options[2] = S3L_SYMMETRIC

The sparse array has symmetric structure and symmetric values.

options[2] = S3L_SYM_STRUCT

The sparse array has symmetric structure but asymmetric values.

options[3] = S3L_NO_PIVOT

Do not use pivoting.

options[3] = S3L_DO_PIVOT

Use pivoting.

roptions is not currently used. It may be used in the future for specifying such parameters as a drop tolerance for pivoting, a threshold value for determining when a block is considered dense, and an amalgamation constant.

setup is an integer that is returned by S3L_sparse_solve upon completion. It describes the sparse linear solution resulting from this call.

Note - If the internal setup created by this call will be used for additional solutions of the linear system, subsequent calls to S3L_sparse_solve would use this setup value. It will also be used in a subsequent call to S3L_sparse_solve_free to free the internal data associated with this sparse linear system solution.

If the call is made from a Fortran program, error status will be in ier.

For detailed descriptions of the Fortran and C bindings for this routine, see the S3L_sparse_solve(3) man page or the corresponding description in the Sun S3L Software Reference Manual.

Examples showing S3L_sparse_solve in use can be found in:

/opt/SUNWhpc/examples/s3l/spsolve/ex_sparse_solve1.c

/opt/SUNWhpc/examples/s3l/spsolve-f/ex_sp_solve1.f

Solving Sparse Linear Systems by an Iterative Method

S3L_gen_iter_solve solves a linear system of equations of the form A*x = y, where A is a sparse Sun S3L array and x and y are both single- or double-precision dense parallel arrays.

Given a general, square sparse matrix A and a right-hand side vector b, S3L_gen_iter_solve solves the linear system of equations Ax = b, using an iterative algorithm, with or without preconditioning.

The first three arguments to S3L_gen_iter_solve are Sun S3L internal array handles that describe the global general sparse matrix A, the rank 1 global array b, and the rank 1 global array x.

The sparse matrix A is produced by a prior call to one of the following sparse routines:

S3L_declare_sparse

S3L_read_sparse

S3L_rand_sparse

S3L_convert_sparse

The global rank 1 arrays, b and x, have the same data type and precision as the sparse matrix A, and both have a length equal to the order of A.

Two local rank 1 arrays, iparm (integer array) and rparm (real array), provide user control over various aspects of S3L_gen_iter_solve behavior, including:

Choice of algorithm to be used

Type of preconditioner to use on A

Flags to select the initial guess to the solution

Maximum number of iterations to be taken by the solver

If the restarted GMRES algorithm is chosen, selection of the size of the Krylov subspace

Tolerance values to be used by the stopping criterion

If the Richardson algorithm is chosen, selection of the scaling factor to be used

The options supported by the iparm and rparm arguments are described in the following subsections.

Note - iparm and rparm must be preallocated and initialized before S3L_gen_iter_solve is called. To enable the default condition for any parameter, set it to 0. Otherwise, initialize the arguments with the appropriate parameter values, as described in the following subsections.

Algorithm

S3L_gen_iter_solve attempts to solve Ax = b using one of the following iterative solution algorithms. The choice of algorithm is determined by the value supplied for the parameter iparm[S3L_iter_solver]. The various options available for this parameter are listed and described in TABLE 13-1.


Option	Description
`S3L_bcgs`	BiConjugate Gradient Stabilized (Bi-CGSTAB)
`S3L_cgs`	Conjugate Gradient Squared (CGS)
`S3L_cg`	Conjugate Gradient (CG)
`S3L_cr`	Conjugate Residuals (CR)
`S3L_gmres`	Generalized Minimum Residual (GMRES) - default
`S3L_qmr`	Quasi-Minimal Residual (QMR)
`S3L_richardson`	Richardson method

Preconditioning

S3L_gen_iter_solve implements left preconditioning. That is, preconditioning is applied to the linear system Ax = b by:

Q^-1 A = Q^-1 b

where Q is the preconditioner and Q^-1 denotes the inverse of Q. The supported preconditioners are listed in TABLE 13-2.


Option	Description
`S3L_none`	No preconditioning will be done (default).
`S3L_jacobi`	Point Jacobi preconditioner will be used. Note that this option is not supported when the sparse matrix `A` is represented under `S3L_SPARSE_VBR` format.
`S3L_bjacobi`	Block Jacobi preconditioner will be used. Note that this option is supported only when the sparse matrix `A` is represented under `S3L_SPARSE_VBR` format.
`S3L_ilu`	Use a simplified ILU(0)--the Incomplete LU factorization of level zero preconditioner. This preconditioner modifies only diagonal nonzero elements of the matrix. Note that this option is not supported when the sparse matrix `A` is represented under `S3L_SPARSE_VBR` format.

Convergence/Divergence Criteria

The iparm[S3L_iter_conv] parameter selects the criterion to be used for stopping computation. Currently, the single valid option for this parameter is S3L_r0, which selects the default criterion for both convergence and divergence. The convergence criterion is satisfied when:

err = ||rj||_2 / ||r0||_2 < epsilon

and the divergence criterion is met when:

err = ||rj||_2 / ||r0||_2 > 10000.0

where:

rj and r0 are the residuals obtained at iterations j and 0.

||.||_2 is the 2-norm.

epsilon is the desired convergence tolerance stored in rparm[S3L_iter_tol].

10000.0 is the divergence tolerance, which is set internally in the solver.

Initial Guess

The parameter iparm[S3L_iter_init] determines the contents of the initial guess to the solution of the linear system as follows:

0 - Applies zero as the initial guess. This is the default.

1 - Applies the value contained in array x as the initial guess. For this case, the user must initialize x before calling S3L_gen_iter_solve.

Maximum Iterations

On input, the iparm[S3L_iter_maxiter] parameter specifies the maximum number of iterations to be taken by the solver. Set to 0 to select the default, which is 10000.

On output, iparm[S3L_iter_maxiter] contains the total number of iterations taken by the solver at the time of termination.

Krylov Subspace

If the restarted GMRES algorithm is selected, iparm[S3L_iter_kspace] specifies the size of the Krylov subspace to be used. The default is 30.

Stopping Criterion Tolerance

On input, rparm[S3L_iter_tol] specifies the tolerance values to be used by the stopping criterion. Its default is 10-8.

On output, rparm[S3L_iter_tol] contains the computed error, err, according to the convergence criteria. See the iparm[S3L_iter_conv] description for details.

Richardson Scaling Factor

If the Richardson method is selected, rparm[S3L_rich_scale] specifies the scaling factor to be used. The default value is 1.0.

Iteration Termination

S3L_gen_iter_solve terminates the iteration when one of the following conditions is met:

The computation has satisfied the convergence criterion.

The computation has diverged.

An algorithmic breakdown has occurred.

The number of iterations has exceeded the supplied value.

S3L_gen_iter_solve has the following argument syntax:

S3L_gen_iter_solve(A, b, x, iparm, rparm, ier)

A is a Sun S3L array handle that describes a global general sparse matrix.

b is a Sun S3L array handle that describes a global array of rank 1. It has the same data type and precision as A. b contains the right-hand side vector of the linear problem.

x is a Sun S3L array handle that describes a global array of rank 1. It has the same data type and precision as A and b. At the start, x contains the initial guess for the solution to the linear system. Upon exit, x contains the converged solution. If the computation breaks down or diverges, x will contain the results of the most recent iteration.

iparm is an integer local array of rank 1 and length s3l_iter_iparm_size. At the start, iparm options have the following uses:

iparm[S3l_iter_solver] - Specifies the iterative algorithm to be used. Set it to 0 to use the default solver GMRES. See the Description section for details.

iparm[S3l_iter_pc] - Specifies the preconditioner to be used. Set it to 0 to use the default option, S3L_none.

iparm[S3l_iter_conv] - Selects the criterion to be used for stopping the computation.

iparm[S3l_iter_init] - Specifies the contents of the initial guess to the solution of the linear system.

iparm[S3l_iter_maxiter] - Specifies the maximum number of iterations to be taken by the solver.

iparm[S3l_iter_kspace] - Specifies the size of the Krylov subspace for restarted GMRES.

Upon exit, iparm contains the total number of iterations taken by the solver at the time of termination.

rparm is a real local array with the same precision as x and a length equal to S3L_iter_rparm_size. At the start, it provides the following options:

rparm[S3L_iter_tol] - Specifies the tolerance values to be used by the stopping criterion. It has a default of 10-8.

rparm[S3L_rich_scale] - Specifies the scaling factor to be used in the Richardson method. The default is 1.0.

Upon exit, rparm contains the value of err, according to the stopping criterion, as computed for the last iteration.

If the call is made from a Fortran program, error status will be in ier.

For detailed descriptions of the Fortran and C bindings for this routine, see the S3L_gen_iter_solve(3) man page or the corresponding description in the Sun S3L Software Reference Manual.

Examples showing S3L_gen_iter_solve in use can be found in:

/opt/SUNWhpc/examples/s3l/iter/ex_iter.c

/opt/SUNWhpc/examples/s3l/iter-f/ex_iter.f

Deallocating a Sparse Linear System Solver

S3L_sparse_solve_free frees all internal data associated with the solution of a sparse linear system.

S3L_sparse_solve_free has the following argument syntax:

S3L_sparse_solve_free(setup, ier)

setup is an integer associated with a particular sparse linear solution. It was previously returned by the S3L_sparse_solve call that produced the solution.

If the call is made from a Fortran program, error status will be in ier.

For detailed descriptions of the Fortran and C bindings for this routine, see the S3L_sparse_solve_free(3) man page or the corresponding description in the Sun S3L Software Reference Manual.

Examples showing S3L_sparse_solve_free in use can be found in:

/opt/SUNWhpc/examples/s3l/spsolve/ex_sparse_solve1.c

/opt/SUNWhpc/examples/s3l/spsolve-f/ex_sp_solvef1.f