Engineering and Scientific Subroutine Library for AIX Version 3 Release 3: Guide and Reference

DSMGCG--General Sparse Matrix Iterative Solve Using Compressed-Matrix Storage Mode

This subroutine solves a general sparse linear system of equations using an iterative algorithm, conjugate gradient squared or generalized minimum residual, with or without preconditioning by an incomplete LU factorization. The subroutine is suitable for positive real matrices--that is, when the symmetric part of the matrix, (A+A^T)/2, is positive definite. The sparse matrix is stored in compressed-matrix storage mode. Matrix A and vectors x and b are used:

Ax = b

where A, x, and b contain long-precision real numbers.

Notes:

These subroutines are provided only for migration purposes. You get better performance and a wider choice of algorithms if you use the DSRIS subroutine.
If your sparse matrix is stored by rows, as defined in Storage-by-Rows, you should first use the utility subroutine DSRSM to convert your sparse matrix to compressed-matrix storage mode. See DSRSM--Convert a Sparse Matrix from Storage-by-Rows to Compressed-Matrix Storage Mode.

Syntax

Fortran	CALL DSMGCG (`m`, `nz`, `ac`, `ka`, `lda`, `b`, `x`, `iparm`, `rparm`, `aux1`, `naux1`, `aux2`, `naux2`)
C and C++	dsmgcg (`m`, `nz`, `ac`, `ka`, `lda`, `b`, `x`, `iparm`, `rparm`, `aux1`, `naux1`, `aux2`, `naux2`);
PL/I	CALL DSMGCG (`m`, `nz`, `ac`, `ka`, `lda`, `b`, `x`, `iparm`, `rparm`, `aux1`, `naux1`, `aux2`, `naux2`);

On Entry

m

is the order of the linear system Ax = b and the number of rows in sparse matrix A. Specified as: a fullword integer; m >= 0.

nz

is the maximum number of nonzero elements in each row of sparse matrix A. Specified as: a fullword integer; nz >= 0.

ac

is the array, referred to as AC, containing the values of the nonzero elements of the sparse matrix, stored in compressed-matrix storage mode. Specified as: an lda by (at least) nz array, containing long-precision real numbers.

ka

is the array, referred to as KA, containing the column numbers of the matrix A elements stored in the corresponding positions in array AC. Specified as: an lda by (at least) nz array, containing fullword integers, where 1 <= (elements of KA) <= m.

lda

is the leading dimension of the arrays specified for ac and ka. Specified as: a fullword integer; lda > 0 and lda >= m.

b

is the vector b of length m, containing the right-hand side of the matrix problem. Specified as: a one-dimensional array of (at least) length m, containing long-precision real numbers.

x

is the vector x of length m, containing your initial guess of the solution of the linear system. Specified as: a one-dimensional array of (at least) length m, containing long-precision real numbers. The elements can have any value, and if no guess is available, the value can be zero.

iparm

is an array of parameters, IPARM(i), where:

IPARM(1) controls the number of iterations.
If IPARM(1) > 0, IPARM(1) is the maximum number of iterations allowed.
If IPARM(1) = 0, the following default values are used:

IPARM(1) = 300
IPARM(2) = 0
IPARM(3) = 10
RPARM(1) = 10^-6
IPARM(2) is the flag used to select the iterative procedure used in this subroutine.
If IPARM(2) = 0, the conjugate gradient squared method is used.
If IPARM(2) = k, the generalized minimum residual method, restarted after k steps, is used. Note that the size of the work area aux1 becomes larger as k increases. A value for k in the range of 5 to 10 is suitable for most problems.
IPARM(3) is the flag that determines whether the system is to be preconditioned by an incomplete LU factorization with no fill-in.
If IPARM(3) = 0, the system is not preconditioned.
If IPARM(3) = 10, the system is preconditioned by an incomplete LU factorization.
If IPARM(3) = -10, the system is preconditioned by an incomplete LU factorization, where the factorization matrix was computed in an earlier call to this subroutine and is stored in aux2.
IPARM(4), see On Return.

Specified as: an array of (at least) length 4, containing fullword integers, where:

IPARM(1) >= 0

IPARM(2) >= 0

IPARM(3) = 0, 10, or -10

rparm

is an array of parameters, RPARM(i), where:

RPARM(1) > 0, is the relative accuracy epsilon used in the stopping criterion. The iterative procedure is stopped when:

||b-Ax||₂ / ||x||₂ < epsilon

RPARM(2) is reserved.

RPARM(3), see On Return.

Specified as: a one-dimensional array of (at least) length 3, containing long-precision real numbers.

aux1

has the following meaning:

If naux1 = 0 and error 2015 is unrecoverable, aux1 is ignored.

Otherwise, it is a storage work area used by this subroutine, which is available for use by the calling program between calls to this subroutine. Its size is specified by naux1.

Specified as: an area of storage, containing long-precision real numbers.

naux1

is the size of the work area specified by aux1--that is, the number of elements in aux1. Specified as: a fullword integer, where:

If naux1 = 0 and error 2015 is unrecoverable, DSMGCG dynamically allocates the work area used by this subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, it must have at least the following value, where:

If IPARM(2) = 0, use naux1 >= 7m.

If IPARM(2) > 0, use naux1 >= (k+2)m+k(k+4)+1, where k = IPARM(2).

aux2

is the storage work area used by this subroutine. If IPARM(3) = -10, aux2 must contain the incomplete LU factorization of matrix A, computed in an earlier call to DSMGCG. The size of aux2 is specified by naux2.

Specified as: an area of storage, containing long-precision real numbers.

naux2

is the size of the work area specified by aux2--that is, the number of elements in aux2. Specified as: a fullword integer. When IPARM(3) = 10, naux2 must have at least the following value: naux2 >= 3+2m+1.5nz(m).

On Return

x

is the vector x of length m, containing the solution of the system Ax = b. Returned as: a one-dimensional array of (at least) length m, containing long-precision real numbers.

iparm

is an array of parameters, IPARM(i), where:

IPARM(1) is unchanged.

IPARM(2) is unchanged.

IPARM(3) is unchanged.

IPARM(4) contains the number of iterations performed by this subroutine.

Returned as: a one-dimensional array of length 4, containing fullword integers.

rparm

is an array of parameters, RPARM(i), where:

RPARM(1) is unchanged.

RPARM(2) is reserved.

RPARM(3) contains the estimate of the error of the solution. If the process converged, RPARM(3) <= RPARM(1)

Returned as: a one-dimensional array of length 3, containing long-precision real numbers.

aux2

is the storage work area used by this subroutine.

If IPARM(3) = 10, aux2 contains the incomplete LU factorization of matrix A.

If IPARM(3) = -10, aux2 is unchanged.

See Notes for additional information on aux2. Returned as: an area of storage, containing long-precision real numbers.

Notes

When IPARM(3) = -10, this subroutine uses the incomplete LU factorization in aux2, computed in an earlier call to this subroutine. When IPARM(3) = 10, this subroutine computes the incomplete LU factorization and stores it in aux2.
If you solve the same sparse linear system of equations several times with different right-hand sides using the preconditioned algorithm, specify IPARM(2) = 10 on the first invocation. The incomplete factorization is stored in aux2. You may save computing time on subsequent calls by setting IPARM(3) equal to -10. In this way, the algorithm reutilizes the incomplete factorization that was computed the first time. Therefore, you should not modify the contents of aux2 between calls.
Matrix A must have no common elements with vectors x and b; otherwise, results are unpredictable.
In the iterative solvers for sparse matrices, the relative accuracy epsilon (RPARM(1)) must be specified "reasonably" (10^-4 to 10^-8). The algorithm computes a sequence of approximate solution vectors x that converge to the solution. The iterative procedure is stopped when the norm of the residual is sufficiently small--that is, when:

||b-Ax||₂ / ||x||₂ < epsilon

As a result, if you specify a larger epsilon, the algorithm takes fewer iterations to converge to a solution. If you specify a smaller epsilon, the algorithm requires more iterations and computer time, but converges to a more precise solution. If the value you specify is unreasonably small, the algorithm may fail to converge within the number of iterations it is allowed to perform.
For a description of how sparse matrices are stored in compressed-matrix storage mode, see Compressed-Matrix Storage Mode.
On output, array AC is not bitwise identical to what it was on input because the matrix A is scaled before starting the iterative process and is unscaled before returning control to the user.
You have the option of having the minimum required value for naux dynamically returned to your program. For details, see Using Auxiliary Storage in ESSL.

Function

The linear system:

Ax = b

is solved using either the conjugate gradient squared method or the generalized minimum residual method, with or without preconditioning by an incomplete LU factorization, where:

A is a sparse matrix of order m, stored in compressed-matrix storage mode in arrays AC and KA.

x is a vector of length m.

b is a vector of length m.

See references [86] and [88]. [36].

If your program uses a sparse matrix stored by rows and you want to use this subroutine, first convert your sparse matrix to compressed-matrix storage mode by using the subroutine DSRSM described on page DSRSM--Convert a Sparse Matrix from Storage-by-Rows to Compressed-Matrix Storage Mode.

Error Conditions

Resource Errors

Error 2015 is unrecoverable, naux1 = 0, and unable to allocate work area.

Computational Errors

The following errors, with their corresponding return codes, can occur in this subroutine. For details on error handling, see What Can You Do about ESSL Computational Errors?.

For error 2110, return code 1 indicates that the subroutine exceeded IPARM(1) iterations without converging. Vector x contains the approximate solution computed at the last iteration.
For error 2111, return code 2 indicates that aux2 contains an incorrect factorization. The subroutine has been called with IPARM(3) = -10, and aux2 contains an incomplete factorization of the input matrix A that was computed by a previous call to the subroutine when IPARM(3) = 10. This error indicates that aux2 has been modified since the last call to the subroutine, or that the input matrix is not the same as the one that was factored. If the default action has been overridden, the subroutine can be called again with the same parameters, with the exception of IPARM(3) = 0 or 10.
For error 2112, return code 3 indicates that the incomplete LU factorization of A could not be completed, because one pivot was 0.
For error 2116, return code 4 indicates that the matrix is singular, because all elements in one row of the matrix contain 0. Array AC is partially modified and does not represent the same matrix as on entry.

Input-Argument Errors

m < 0
lda < 1
lda < m
nz < 0
nz = 0 and m > 0
IPARM(1) < 0
IPARM(2) < 0
IPARM(3) <> 0, 10, or -10
RPARM(1) < 0
RPARM(2) < 0
Error 2015 is recoverable or naux1<>0, and naux1 is too small--that is, less than the minimum required value. Return code 5 is returned if error 2015 is recoverable.
naux2 is too small--that is, less than the minimum required value. Return code 5 is returned if error 2015 is recoverable.

Example 1

This example finds the solution of the linear system Ax = b for the sparse matrix A, which is stored in compressed-matrix storage mode in arrays AC and KA. The system is solved using the conjugate gradient squared method. Matrix A is:

        *                                                   *
        | 2.0  0.0   0.0  0.0   0.0   0.0   0.0   0.0   0.0 |
        | 0.0  2.0  -1.0  0.0   0.0   0.0   0.0   0.0   0.0 |
        | 0.0  1.0   2.0  0.0   0.0   0.0   0.0   0.0   0.0 |
        | 1.0  0.0   0.0  2.0  -1.0   0.0   0.0   0.0   0.0 |
        | 0.0  0.0   0.0  1.0   2.0  -1.0   0.0   0.0   0.0 |
        | 0.0  0.0   0.0  0.0   1.0   2.0  -1.0   0.0   0.0 |
        | 0.0  0.0   0.0  0.0   0.0   1.0   2.0  -1.0   0.0 |
        | 0.0  0.0   0.0  0.0   0.0   0.0   1.0   2.0  -1.0 |
        | 0.0  0.0   0.0  0.0   0.0   0.0   0.0   1.0   2.0 |
        *                                                   *

Note:: For input matrix KA, ( . ) indicates any value between 1 and 9.

Call Statement and Input

             M   NZ   AC   KA   LDA   B   X   IPARM   RPARM   AUX1   NAUX1   AUX2   NAUX2
             |    |    |    |    |    |   |     |       |      |       |      |       |
CALL DSMGCG( 9 ,  3 , AC , KA ,  9  , B , X , IPARM , RPARM , AUX1 ,  63   , AUX2 ,   0 )

IPARM(1) =  20
IPARM(2) =  0
IPARM(3) =  0
RPARM(1) =  1.D-7

        *                 *
        | 2.0   0.0   0.0 |
        | 2.0  -1.0   0.0 |
        | 1.0   2.0   0.0 |
        | 1.0   2.0  -1.0 |
AC   =  | 1.0   2.0  -1.0 |
        | 1.0   2.0  -1.0 |
        | 1.0   2.0  -1.0 |
        | 1.0   2.0  -1.0 |
        | 1.0   2.0   0.0 |
        *                 *

        *         *
        | 1  .  . |
        | 2  3  . |
        | 2  3  . |
        | 1  4  5 |
KA   =  | 4  5  6 |
        | 5  6  7 |
        | 6  7  8 |
        | 7  8  9 |
        | 8  9  . |
        *         *

B        =  (2.0, 1.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0)
X        =  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

Output

X        =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
IPARM(4) =  9
RPARM(3) =  0.150D-19

Example 2

This example finds the solution of the linear system Ax = b for the same sparse matrix A as in Example 1, which is stored in compressed-matrix storage mode in arrays AC and KA. The system is solved using the generalized minimum residual method, restarted after 5 steps and preconditioned with an incomplete LU factorization. Most of the input is the same as in Example 1.

Note:: For input matrix KA, ( . ) indicates any value between 1 and 9.

Call Statement and Input

             M   NZ   AC   KA   LDA   B   X   IPARM   RPARM   AUX1   NAUX1   AUX2   NAUX2
             |    |    |    |    |    |   |     |       |      |       |      |       |
CALL DSMGCG( 9 ,  3 , AC , KA ,  9  , B , X , IPARM , RPARM , AUX1 ,  109  , AUX2 ,  46 )

IPARM(1) = 20
IPARM(2) = 5
IPARM(3) = 10
RPARM(1) = 1.D-7
AC =(same as input AC in Example 1)
KA =(same as input KA in Example 1)
B = (2.0, 1.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0)
X = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

Output

X        =  (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
IPARM(4) =  2
RPARM(3) =  0.290D-15

[ Top of Page | Previous Page | Next Page | Table of Contents | Index ]