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

DGSF--General Sparse Matrix Factorization Using Storage by Indices, Rows, or Columns

This subroutine factors sparse matrix A by Gaussian elimination, using a modified Markowitz count with threshold pivoting. The sparse matrix can be stored by indices, rows, or columns. To solve the system of equations, follow the call to this subroutine with a call to DGSS.

Syntax

On Entry

iopt

indicates the storage technique used for sparse matrix A, where:

If iopt = 0, it is stored by indices.

If iopt = 1, it is stored by rows.

If iopt = 2, it is stored by columns.

Specified as: a fullword integer; iopt = 0, 1, or 2.

n

is the order n of sparse matrix A. Specified as: a fullword integer; n >= 0.

nz

is the number of elements in sparse matrix A, stored in an array, referred to as A. Specified as: a fullword integer; nz > 0.

a

is the sparse matrix A, to be factored, stored in an array, referred to as A. Specified as: an array of length lna, containing long-precision real numbers.

ia

is the array, referred to as IA, where:

If iopt = 0, it contains the row numbers that correspond to the elements in array A.

If iopt = 1, it contains the row pointers.

If iopt = 2, it contains the row numbers that correspond to the elements in array A.

Specified as: an array of length lna, containing fullword integers; IA(i) >= 1. See Sparse Matrix for more information on storage techniques.

ja

is the array, referred to as JA, where:

If iopt = 0, it contains the column numbers that correspond to the elements in array A.

If iopt = 1, it contains the column numbers that correspond to the elements in array A.

If iopt = 2, it contains the column pointers.

Specified as: an array of length lna, containing fullword integers; JA(i) >= 1. See Sparse Matrix for more information on storage techniques.

lna

is the length of the arrays specified for a, ia, and ja. Specified as: a fullword integer; lna > 2nz. If you do not specify a sufficient amount, it results in an error. See Error Conditions.

The size of lna depends on the structure of the input matrix. The requirement that lna > 2nz does not guarantee a successful run of the program. If the input matrix is expected to have many fill-ins, lna should be set larger. Larger lna may result in a performance improvement.

For details on how lna relates to storage compressions, see Performance and Accuracy Considerations.

iparm

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

• IPARM(1) determines whether the default values for iparm and rparm are used by this subroutine.

  If IPARM(1) = 0, the following default values are used:

  IPARM(2) = 10

  IPARM(3) = 1

  IPARM(4) = 0

  RPARM(1) = 10^-12

  RPARM(2) = 0.1

  If IPARM(1) = 1, the default values are not used.

• IPARM(2) determines the number of minimal Markowitz counts that are examined to determine a pivot. (See reference [101].)
• IPARM(3) has the following meaning, where:

  If IPARM(3) = 0, this subroutine checks the values in arrays IA and JA.

  If IPARM(3) = 1, this subroutine assumes that the input values are correct in arrays IA and JA.

• IPARM(4) has the following meaning, where:

  If IPARM(4) = 0, this computation is not performed.

  If IPARM(4) = 1, this subroutine computes:

  The absolute value of the smallest pivot element

  The absolute value of the largest element in U.

  These values are stored in OPARM(2) and OPARM(3), respectively.

• IPARM(5) is reserved.

Specified as: an array of (at least) length 5, containing fullword integers, where the iparm values must be:

IPARM(1) = 0 or 1

IPARM(2) >= 1

IPARM(3) = 0 or 1

IPARM(4) = 0 or 1

rparm

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

• RPARM(1) contains the lower bound of the absolute value of all elements in the matrix. If a pivot element is less than this number, the matrix is reported as singular. Any computed element whose absolute value is less than this number is set to 0.
• RPARM(2) is the threshold pivot tolerance used to control the choice of pivots.
• RPARM(3) is reserved.
• RPARM(4) is reserved.
• RPARM(5) is reserved.

Specified as: a one-dimensional array of (at least) length 5, containing long-precision real numbers, where the rparm values must be:

RPARM(1) >= 0.0

0.0 <= RPARM(2) <= 1.0

For additional information about rparm, see Performance and Accuracy Considerations.

oparm

See On Return.

aux

is the storage work area used by this subroutine. Its size is specified by naux. Specified as: an area of storage, containing long-precision real numbers.

naux

is the size of the work area specified by aux--that is, the number of elements in aux. Specified as: a fullword integer; naux >= 10n+100.

On Return

a

is the transformed array, referred to as A, containing the factored matrix A, required as input to DGSS. Returned as: a one-dimensional array of length lna, containing long-precision real numbers.

ia

is the transformed array, referred to as IA, required as input to DGSS. Returned as: a one-dimensional array of length lna, containing fullword integers.

ja

is the transformed array, referred to as JA, required as input to DGSS. Returned as: a one-dimensional array of length lna, containing fullword integers.

oparm

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

• OPARM(1) is the amount of fill-ins for the sparse processing portion of the algorithm.
• OPARM(2) contains the absolute value of the smallest pivot element of the matrix. This value is computed and set only if IPARM(4) = 1.
• OPARM(3) contains the absolute value of the largest element encountered in U after the factorization. This value is computed and set only if IPARM(4) = 1.
• OPARM(4) is reserved.
• OPARM(5) is reserved.

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

aux

is the storage work area used by this subroutine. It contains the information required as input for DGSS. Specified as: an area of storage, containing long-precision real numbers.

Notes

1. For a description of the three storage techniques used by this subroutine for sparse matrices, see Sparse Matrix.
2. 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 matrix A is factored by Gaussian elimination, using a modified Markowitz count with threshold pivoting to compute the sparse LU factorization of A:

LU = PAQ

where:

A is a general sparse matrix of order n, stored by indices, columns, or rows in arrays A, IA, and JA.

L is a unit lower triangular matrix.

U is an upper triangular matrix.

P is a permutation matrix.

Q is a permutation matrix.

To solve the system of equations, follow the call to this subroutine with a call to DGSS. If n is 0, no computation is performed. See references [10], [47], and [93].

Error Conditions

Computational Errors

1. If this subroutine has to perform storage compressions, an attention message is issued. When this occurs, the performance of this subroutine is affected. The performance can be improved by increasing the value specified for lna.
2. The following errors with their corresponding return codes can occur in this subroutine. Where a value of i is indicated, it can be determined at run time by use of the ESSL error-handling facilities. To obtain this information, you must use ERRSET to change the number of allowable errors for that particular error code in the ESSL error option table; otherwise, the default value causes your program to terminate when the error occurs. For details, see What Can You Do about ESSL Computational Errors?.
   • For error 2117, return code 2 indicates that the pivot element in a column, i, is smaller than the value specified in RPARM(1).
   • For error 2118, return code 3 indicates that pivot element in a row, i, is smaller than the value specified in RPARM(1).
   • For error 2120, return code 4 indicates that a row, i, is found empty on factorization. The matrix is singular.
   • For error 2121, return code 5 indicates that a column is found empty on factorization. The matrix is singular.
   • For error 2119, return code 6 indicates that the storage space indicated by lna is insufficient.
   • For error 2122, return code 7 indicates that no pivot element was found in the active submatrix.

Input-Argument Errors

1. iopt <> 0, 1, or 2
2. n < 0
3. nz <= 0
4. lna <= 2nz
5. IPARM(1) <> 0 or 1
6. IPARM(2) <= 0
7. IPARM(3) <> 0 or 1
8. IPARM(4) <> 0 or 1
9. RPARM(1) < 0.0
10. RPARM(2) < 0.0 or RPARM(2) > 1.0
11. iopt = 1 and ia(i) >= ia (i+1), i = 1, n
12. iopt = 2 and ja(i) >= ja(i+1), i = 1, n
13. iopt = 0 or 1 and ja(i) < 1 or ja(i) > n, i = 1, nz
14. iopt = 0 or 1 and ia(i) < 1 or ia(i) > n, i = 1, nz
15. There are duplicate indices in a row or column of the input matrix.
16. The matrix is singular if a row or column of the input matrix is empty.
17. naux is too small--that is, less than the minimum required value. Return code 1 is returned if error 2015 is recoverable.

Example

This example factors 5 by 5 sparse matrix A, which is stored by indices in arrays A, IA, and JA. The three storage techniques are shown in this example, and the output is the same regardless of the storage technique used. The matrix is factored using Gaussian elimination with threshold pivoting. Matrix A is:

                       *                         *
                       | 2.0  0.0  4.0  0.0  0.0 |
                       | 1.0  1.0  0.0  0.0  3.0 |
                       | 0.0  0.0  3.0  4.0  0.0 |
                       | 2.0  2.0  0.0  1.0  5.0 |
                       | 0.0  0.0  1.0  1.0  0.0 |
                       *                         *

Note:

In this example, only nonzero elements are used as input to the matrix.

Call Statement and Input (Storage-By-Indices)

          IOPT  N  NZ  A  IA  JA  LNA  IPARM  RPARM  OPARM  AUX  NAUX
           |    |   |  |   |   |   |     |      |      |     |    |
CALL DGSF( 0  , 5, 13, A, IA, JA, 27 , IPARM, RPARM, OPARM, AUX, 150 )

A        =  (2.0, 1.0, 1.0, 3.0, 4.0, 1.0, 5.0, 2.0, 2.0, 1.0, 1.0,
             4.0, 3.0, . , . , . , . , . , . , . , . , . , . , . , . ,
             . , . )
IA       =  (1, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 1, 2, . , . , . , . ,
             . , . , . , . , . , . , . , . , . , . )
JA       =  (1, 1, 2, 3, 4, 4, 5, 1, 2, 3, 4, 3, 5, . , . , . , . ,
             . , . , . , . , . , . , . , . , . , . )
IPARM    =  (1, 3, 1, 1)
RPARM    =  (1.D-12, 0.1D0)

Call Statement and Input (Storage-By-Rows)

          IOPT  N  NZ  A  IA  JA  LNA  IPARM  RPARM  OPARM  AUX  NAUX
           |    |   |  |   |   |   |     |      |      |     |    |
CALL DGSF( 1  , 5, 13, A, IA, JA, 27 , IPARM, RPARM, OPARM, AUX, 150 )

A        =  (2.0, 4.0, 1.0, 1.0, 3.0, 3.0, 4.0, 2.0, 2.0, 1.0, 5.0,
             1.0, 1.0, . , . , . , . , . , . , . , . , . , . , . , . ,
             . , . )
IA       =  (1, 3, 6, 8, 12, 14, . , . , . , . , . , . , . , . , . ,
             . , . , . , . , . )
JA       =  (1, 3, 1, 2, 5, 3, 4, 1, 2, 4, 5, 3, 4, . , . , . , . ,
             . , . , . , . , . , . , . , . , . , . )
IPARM    =  (1, 3, 1, 1)
RPARM    =  (1.D-12, 0.1D0)

Call Statement and Input (Storage-By-Columns)

          IOPT  N  NZ  A  IA  JA  LNA  IPARM  RPARM  OPARM  AUX  NAUX
           |    |   |  |   |   |   |     |      |      |     |    |
CALL DGSF( 2  , 5, 13, A, IA, JA, 27 , IPARM, RPARM, OPARM, AUX, 150 )

A        =  (2.0, 1.0, 2.0, 1.0, 2.0, 4.0, 3.0, 1.0, 4.0, 1.0, 1.0,
             3.0, 5.0, . , . , . , . , . , . , . , . , . , . , . , . ,
             . , . )
IA       =  (1, 2, 4, 2, 4, 1, 3, 5, 3, 4, 5, 2, 4, . , . , . , . ,
             . , . , . , . , . , . , . , . , . , . )
JA       =  (1, 4, 6, 9, 12, 14, . , . , . , . , . , . , . , . , . ,
             . , . , . , . , . )
IPARM    =  (1, 3, 0, 1)
RPARM    =  (1.D-12, 0.1D0)

Output

A        =  (0.5, . , 0.3, 1.0, . , 1.0, . , 3.0, . , . , . , 1.0,
             1.0, . , . , . , . , . , . , . , -1.7, -0.5, -1.0, -1.0,
             4.0, -3.0, -4.0)
IA       =  (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, . , . , . , . ,
             . , . , . , 2, 1, 1, 3, 3, 5, 5)
JA       =  (1, 0, 5, 2, 0, 4, 0, 2, 0, 0, 0, 3, 4, . , . , . , . ,
             . , . , . , 4, 2, 4, 4, 1, 3, 1)
OPARM    =  (1.000000, 0.333333, 3.000000)