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

DSDGCG--General Sparse Matrix Iterative Solve Using Compressed-Diagonal 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-diagonal storage mode. Matrix A and vectors x and b are used:

Ax = b

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

Syntax

Fortran	CALL DSDGCG (`m`, `nd`, `ad`, `lda`, `la`, `b`, `x`, `iparm`, `rparm`, `aux1`, `naux1`, `aux2`, `naux2`)
C and C++	dsdgcg (`m`, `nd`, `ad`, `lda`, `la`, `b`, `x`, `iparm`, `rparm`, `aux1`, `naux1`, `aux2`, `naux2`);
PL/I	CALL DSDGCG (`m`, `nd`, `ad`, `lda`, `la`, `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.

nd

is the number of nonzero diagonals stored in the columns of array AD, the number of columns in array AD, and the number of elements in array LA. Specified as: a fullword integer; it must have the following value, where:

If m > 0, then nd > 0.

If m = 0, then nd >= 0.

ad

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

lda

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

la

is the array, referred to as LA, containing the diagonal numbers k for the diagonals stored in each corresponding column in array AD. For an explanation of how diagonal numbers are stored, see Compressed-Diagonal Storage Mode.

Specified as: a one-dimensional array of (at least) length nd, containing fullword integers, where 1-m <= (elements of LA) <= (m-1).

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:

If 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, DSDGCG dynamically allocates the work area used by this subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, naux1 > 0 and 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 = PARM(2).

aux2

is a 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 DSDGCG. 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 or -10, naux2 must have at least the following value: naux2 >= 3+2m+1.5nd(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(3) = 10 on the first invocation. The incomplete factorization is stored in aux2. You may save computing time on subsequent calls by setting IPARM(3) = -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-diagonal storage mode, see Compressed-Diagonal Storage Mode.
On output, array AD is not bitwise identical to what it was on input, because 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-diagonal storage mode in arrays AD and LA.

x is a vector of length m.

b is a vector of length m.

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

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
nd < 0
nd = 0 and m > 0
IPARM(1) < 0
IPARM(2) < 0
IPARM(3) <> 0, 10, or -10
RPARM(1) < 1.D0
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-diagonal storage mode in arrays AD and LA. The system is solved using the conjugate gradient squared method. Matrix A is:

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

Call Statement and Input

             M   ND   AD   LDA   LA   B   X   IPARM   RPARM   AUX1   NAUX1   AUX2   NAUX2
             |    |    |    |     |   |   |     |       |      |       |      |       |
CALL DSDGCG( 9 ,  3 , AD ,  9  , LA , B , X , IPARM , RPARM , AUX1 ,  63   , AUX2 ,   0 )

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

        *                *
        | 2.0  -1.0  0.0 |
        | 2.0  -1.0  0.0 |
        | 2.0  -1.0  0.0 |
        | 2.0  -1.0  0.0 |
AD   =  | 2.0  -1.0  1.0 |
        | 2.0  -1.0  1.0 |
        | 2.0  -1.0  1.0 |
        | 2.0   0.0  1.0 |
        | 2.0   0.0  1.0 |
        *                *

LA       =  (0, 2, -4)
B        =  (1, 1, 1, 1, 2, 2, 2, 3, 3)
X        =  (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) =  8
RPARM(3) =  0.308D-17

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-diagonal storage mode in arrays AD and LA. 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.

Call Statement and Input

             M   ND   AD   LDA   LA   B   X   IPARM   RPARM   AUX1   NAUX1   AUX2   NAUX2
             |    |    |    |     |   |   |     |       |      |       |      |       |
CALL DSDGCG( 9 ,  3 , AD ,  9  , LA , B , X , IPARM , RPARM , AUX1 ,  109  , AUX2 ,  46 )

IPARM(1) = 20
IPARM(2) = 5
IPARM(3) = 10
RPARM(1) = 1.D-7
AD =(same as input AD in Example 1)
LA =(same as input LA in Example 1)
B = (1, 1, 1, 1, 2, 2, 2, 3, 3)
X = (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) =  6
RPARM(3) =  0.250D-15

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