Iterative Solver

S3L_gen_iter_solve

Description

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 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

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 and rparm, 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 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.

iparm is an integer array and rparm is a real array. The options supported by these arguments are described in the subsections titled: "Algorithm," "Preconditioning," "Initial Guess," "Maximum Iterations," "Krylov Subspace," "Stopping Criterion Tolerance," and "Richardson Scaling Factor." The "Iteration Termination" subsection identifies the conditions under which S3L_gen_iter_solve will terminate anoperation.

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 them 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 8-12

Preconditioning

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

Example 8-81

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 8-13.

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.

Syntax

The C and Fortran syntax for S3L_gen_iter_solve are shown below.

C/C++ Syntax

Example 8-82

#include <s3l/s3l-c.h>
#include <s3l/s3l_errno-c.h>
int
S3L_gen_iter_solve(A, b, x, iparm, rparm)
    S3L_array_t        A
    S3L_array_t        b
    S3L_array_t        x
    int                *iparm
    <type>            *rparm

F77/F90 Syntax

Example 8-83

include `s3l/s3l-f.h'
include `s3l/s3l_errno-f.h'
subroutine
S3L_gen_iter_solve(A, b, x, iparm, rparm, ier)
    integer*8          A
    integer*8          b
    integer*8          x
    integer*4         iparm(*)
    <type>            rparm(*)
    integer*4         ier

where <type> is real*4 or real*8 for both C/C++ and F77/F90.

Input

• A - S3L internal array handle for the global general sparse matrix. It is produced by a prior call to one of the following sparse routines:

  • S3L_declare_sparse

  • S3L_read_sparse

  • S3L_rand_sparse

• b - Global array of rank 1, with the same data type and precision as A and x and a length equal to the order of the sparse matrix. b contains the right-hand side vector of the linear problem.

• x - Global array of rank 1, with the same data type and precision as A and b and a length equal to the order of the sparse matrix. On input, x contains the initial guess for the solution to the linear system. Upon successful completion, x contains the converged solution (see the Output section).

• iparm - Integer local array of rank 1 and length s3l_iter_iparm_size, where:

  • iparm[S3l_iter_solver] - Specifies the iterative algorithm to be used. Set it to 0 to use the default solver GMRES. See the Desctription sectino 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.

• rparm - Specifies options for computing all or part of the matrix U.

Output

This function uses the following arguments for output:

• x - Upon successful completion, x contains the converged solution. If the computation breaks down or diverges, x will contain the solution produced by the most recent iteration.

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

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

• ier (Fortran only) - When called from a Fortran program, this function returns error status in ier.

Error Handling

On success, S3L_gen_iter_solve returns S3L_SUCCESS.

S3L_gen_iter_solve performs generic checking of the arrays it accepts as arguments. If an array argument contains an invalid or corrupted value, the function terminates and an error code indicating which value of the array handle was invalid is returned. See Appendix A of this manual for a detailed list of these error codes.

On error, it returns one of the following codes, which are organized by error type.

Input Errors

• S3L_ERR_ARG_NULL - Invalid array x or b or sparse matrix A. They all must be preallocated S3L arrays or sparse matrix.

• S3L_ERR_ARRNOTSQ - Invalid matrix size. Matrix A must be square.

• S3L_ERR_ARG_RANK - Invalid rank for arrays x and b. Both must be rank 1 arrays.

• S3L_ERR_MATCH_DTYPE - x, b, and A do not have the same data type.

• S3L_ERR_MATCH_EXTENTS - The lengths of x and b do not match the size of sparse matrix A. Both must be equal to the order of A.

• S3L_ERR_PARAM_INVALID - Invalid input for iparm or rparm. Both must be preallocated and initialized with the predefined values described in the Description section or set to 0 for the default value.

Computational Errors

• S3L_ERR_ILU_ZRPVT - Encountered a zero pivot in the course of ILU preconditioning.

• S3L_ERR_JACOBI_ZRDIAG - Encountered a zero diagonal in the course of Jacobi preconditioning.

• S3L_ERR_DIVERGE - Computation has diverged.

• S3L_ERR_ITER_BRKDWN - A breakdown has occurred.

• S3L_ERR_MAXITER - The number of iterations has exceeded the value supplied in iparm[S3L_iter_maxiter].

Examples

../examples/s3l/iter/ex_iter.c
../examples/s3l/iter-f/ex_iter.f

Related Functions

S3L_declare_sparse(3)
S3L_read_sparse(3)
S3L_rand_sparse(3)