Integer, Intent (In)	::	n, irowix(*), nzlmx, nzumx
Integer, Intent (Inout)	::	iprm(7*n), nzlumx, ifail
Integer, Intent (Out)	::	il(7n+nzlmx+4), iu(2n+nzumx+1), nnzl, nnzu
Real (Kind=nag_wp), Intent (In)	::	a(*), thresh
Real (Kind=nag_wp), Intent (Out)	::	lval(nzlumx), uval(nzumx), flop

C Header Interface

#include <nag.h>

void	f11mef_ (const Integer n, const Integer irowix[], const double a[], Integer iprm[], const double thresh, const Integer nzlmx, Integer nzlumx, const Integer nzumx, Integer il[], double lval[], Integer iu[], double uval[], Integer nnzl, Integer nnzu, double flop, Integer *ifail)

The routine may be called by the names f11mef or nagf_sparse_direct_real_gen_lu.

3 Description

Given a real sparse matrix

A

, f11mef computes an

L U

factorization of

A

with partial pivoting,

P_{r} A P_{c} = L U

, where

P_{r}

is a row permutation matrix (computed by f11mef),

P_{c}

is a (supplied) column permutation matrix,

L

is unit lower triangular and

U

is upper triangular. The column permutation matrix,

P_{c}

, 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

L U

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 arguments (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.

4 References

Demmel J W, Eisenstat S C, Gilbert J R, Li X S and Li J W H (1999) A supernodal approach to sparse partial pivoting SIAM J. Matrix Anal. Appl. 20 720–755

Demmel J W, Gilbert J R and Li X S (1999) An asynchronous parallel supernodal algorithm for sparse gaussian elimination SIAM J. Matrix Anal. Appl. 20 915–952

5 Arguments

1: $n$ – Integer Input: On entry: $n$ , the order of the matrix $A$ .

Constraint: $n \geq 0$ .
2: $irowix (*)$ – Integer array Input: 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$ . See Section 2.1.3 in the F11 Chapter Introduction.
3: $a (*)$ – Real (Kind=nag_wp) array Input: 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 \times n)$ – Integer array Input/Output: On entry: contains the column permutation which defines the permutation $P_{c}$ and associated data structures as computed by routine f11mdf.

On exit: part of the array is modified to record the row permutation $P_{r}$ 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 $|A_{j j}| \geq t (\max_{i \geq j} |A_{i j}|)$ , use $A_{j j}$ as a pivot, otherwise use $\max_{i \geq j} |A_{i j}|$ . 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 \leq thresh \leq 1.0$ .
6: $nzlmx$ – Integer Input: On entry: indicates the available size of array il. The dimension of il should be at least $7 \times n + nzlmx + 4$ . A good range for nzlmx that works for many problems is $nnz$ to $8 \times 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 \geq 1$ .
7: $nzlumx$ – Integer Input/Output: On entry: indicates the available size of array lval. The dimension of lval should be at least nzlumx.

Constraint: $nzlumx \geq 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$ – Integer Input: On entry: indicates the available sizes of arrays iu and uval. The dimension of iu should be at least $2 \times 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 \times 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.

Constraint: $nzumx \geq 1$ .
9: $il (7 \times n + nzlmx + 4)$ – Integer array Output: On exit: encapsulates the sparsity pattern of matrix $L$ .
10: $lval (nzlumx)$ – Real (Kind=nag_wp) array Output: On exit: records the nonzero values of matrix $L$ and some of the nonzero values of matrix $U$ .
11: $iu (2 \times n + nzumx + 1)$ – Integer array Output: On exit: encapsulates the sparsity pattern of matrix $U$ .
12: $uval (nzumx)$ – Real (Kind=nag_wp) array Output: On exit: records some of the nonzero values of matrix $U$ .
13: $nnzl$ – Integer Output: On exit: the number of nonzero values in the matrix $L$ .
14: $nnzu$ – Integer Output: 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$ – Integer Input/Output: On entry: ifail must be set to $0$ , $- 1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $- 1$ means that an error message is printed while a value of $1$ means that it is not.

If halting is not appropriate, the value $- 1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $0$ is recommended. 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

- 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 = 〈value〉$ .
Constraint: $n \geq 0$ .

On entry, $nzlmx = 〈value〉$ .
Constraint: $nzlmx \geq 1$ .

On entry, $nzlumx = 〈value〉$ .
Constraint: $nzlumx \geq 1$ .

On entry, $nzumx = 〈value〉$ .
Constraint: $nzumx \geq 1$ .

On entry, $thresh = 〈value〉$ .
Constraint: $0.0 \leq thresh \leq 1.0$ .

$ifail = 2$: Insufficient nzlmx.

$ifail = 3$: Insufficient nzumx.

$ifail = 4$: Insufficient nzlumx.

$ifail = 5$: The matrix is singular – no factorization possible.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The computed factors

L

and

U

are the exact factors of a perturbed matrix

A + E

, where

|E| \leq c (n) ε |L| |U|,

c (n)

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 Parallelism and Performance

f11mef is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

f11mef makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

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

A

A call to f11mef may be followed by calls to the routines:

f11mff to solve $A X = B$ or $A^{T} X = B$ ;
f11mgf to estimate the condition number of $A$ ;
f11mmf to estimate the reciprocal pivot growth of the $L U$ factorization.

10 Example

This example computes the

L U

factorization of the matrix

A

, where

A = (\begin{matrix} 2.00 & 1.00 & 0 & 0 & 0 \\ 0 & 0 & 1.00 & - 1.00 & 0 \\ 4.00 & 0 & 1.00 & 0 & 1.00 \\ 0 & 0 & 0 & 1.00 & 2.00 \\ 0 & - 2.00 & 0 & 0 & 3.00 \end{matrix}) .

f11me: FL CL CPP AD

NAG FL Interfacef11mef (direct_​real_​gen_​lu)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
f11mef (direct_real_gen_lu)