3 Description

Given a real sparse matrix A, F11MEF computes an LU factorization of A with partial pivoting, PrAPc=LU, where Pr is a row permutation matrix (computed by F11MEF), Pc is a (supplied) column permutation matrix, L is unit lower triangular and U is upper triangular. The column permutation matrix, Pc, must be computed by a prior call to F11MDF. The matrix A must be presented in the column permuted, compressed column (Harwell–Boeing) format.

The LU factorization is output in the form of four one-dimensional arrays: integer arrays IL and IU and real-valued arrays LVAL and UVAL. These describe the sparsity pattern and numerical values in the L and U matrices. The minimum required dimensions of these arrays cannot be given as a simple function of the size parameters (order and number of nonzero values) of the matrix A. This is due to unpredictable fill-in created by partial pivoting. F11MEF will, on return, indicate which dimensions of these arrays were not adequate for the computation or (in the case of one of them) give a firm bound. You should then allocate more storage and try again.

5 Parameters

Note: the dimension of the array IROWIX
must be at least
nnz, the number of nonzeros of the sparse matrix A.

On entry: the row index array of sparse matrix A.

3: A(*) – REAL (KIND=nag_wp) arrayInput

Note: the dimension of the array A
must be at least
nnz, the number of nonzeros of the sparse matrix A.

On entry: the array of nonzero values in the sparse matrix A.

4: IPRM(7×N) – INTEGER arrayInput/Output

On entry: contains the column permutation which defines the permutation Pc and associated data structures as computed by routine F11MDF.

On exit: part of the array is modified to record the row permutation Pr determined by pivoting.

5: THRESH – REAL (KIND=nag_wp)Input

On entry: the diagonal pivoting threshold, t. At step j of the Gaussian elimination, if Ajj≥t⁢maxi≥jAij, use Ajj as a pivot, otherwise use maxi≥jAij. A value of t=1 corresponds to partial pivoting, a value of t=0 corresponds to always choosing the pivot on the diagonal (unless it is zero).

Suggested value:
THRESH=1.0. Smaller values may result in a faster factorization, but the benefits are likely to be small in most cases. It might be possible to use THRESH=0.0 if you are confident about the stability of the factorization, for example, if A is diagonally dominant.

Constraint:
0.0≤THRESH≤1.0.

6: NZLMX – INTEGERInput

On entry: indicates the available size of array IL. The dimension of IL should be at least 7×N+NZLMX+4. A good range for NZLMX that works for many problems is nnz to 8×nnz, where nnz is the number of nonzeros in the sparse matrix A. If, on exit, IFAIL=2, the given NZLMX was too small and you should attempt to provide more storage and call the routine again.

Constraint:
NZLMX≥1.

7: NZLUMX – INTEGERInput/Output

On entry: indicates the available size of array LVAL. The dimension of LVAL should be at least NZLUMX.

Constraint:
NZLUMX≥1.

On exit: if IFAIL=4, the given NZLUMX was too small and is reset to a value that will be sufficient. You should then provide the indicated storage and call the routine again.

8: NZUMX – INTEGERInput

On entry: indicates the available sizes of arrays IU and UVAL. The dimension of IU should be at least 2×N+NZUMX+1 and the dimension of UVAL should be at least NZUMX. A good range for NZUMX that works for many problems is nnz to 8×nnz, where nnz is the number of nonzeros in the sparse matrix A. If, on exit, IFAIL=3, the given NZUMX was too small and you should attempt to provide more storage and call the routine again.

On exit: records the nonzero values of matrix L and some of the nonzero values of matrix U.

11: IU(2×N+NZUMX+1) – INTEGER arrayOutput

On exit: encapsulates the sparsity pattern of matrix U.

12: UVAL(NZUMX) – REAL (KIND=nag_wp) arrayOutput

On exit: records some of the nonzero values of matrix U.

13: NNZL – INTEGEROutput

On exit: the number of nonzero values in the matrix L.

14: NNZU – INTEGEROutput

On exit: the number of nonzero values in the matrix U.

15: FLOP – REAL (KIND=nag_wp)Output

On exit: the number of floating point operations performed.

16: IFAIL – INTEGERInput/Output

On entry: IFAIL must be set to 0, -1​ or ​1. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value -1​ or ​1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is 0. When the value -1​ or ​1 is used it is essential to test the value of IFAIL on exit.

On exit: IFAIL=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry IFAIL=0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF).

Errors or warnings detected by the routine:

IFAIL=1

On entry,

N<0,

or

NZLMX<1,

or

NZLUMX<1,

or

NZUMX<1,

or

THRESH<0.0,

or

THRESH>1.0.

IFAIL=2

NZLMX was not large enough. You should repeat the call with a larger value of NZLMX, providing more storage for the output array IL.

IFAIL=3

NZUMX was not large enough. You should repeat the call with a larger value of NZUMX, providing more storage for the output arrays IU and UVAL.

IFAIL=4

NZLUMX was not large enough. You should repeat the call with the value of NZLUMX returned on exit, providing more storage for the output array LVAL.

IFAIL=5

The matrix A is singular and no factorization will be attempted.

IFAIL=301

Unable to allocate required internal workspace.

7 Accuracy

The computed factors L and U are the exact factors of a perturbed matrix A+E, where

E≤cnεLU,

cn is a modest linear function of n, and ε is the machine precision, when partial pivoting is used. If no partial pivoting is used, the factorization accuracy can be considerably worse. A call to F11MMF after F11MEF can help determine the quality of the factorization.

8 Further Comments

The total number of floating point operations depends on the sparsity pattern of the matrix A.