NAG Library Routine Document

f01brf  (real_gen_sparse_lu)

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

1:     $\mathbf{n}$ – IntegerInput

On entry: $n$, the order of the matrix $A$.

Constraint: $\mathbf{n}>0$.

2:     $\mathbf{nz}$ – IntegerInput

On entry: the number of nonzero elements in the matrix $A$.

Constraint: $\mathbf{nz}>0$.

3:     $\mathbf{a}\left(\mathbf{licn}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: $\mathbf{a}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{nz}$, must contain the nonzero elements of the sparse matrix $A$. They can be in any order since f01brf will reorder them.

On exit: the nonzero elements in the $LU$ factorization. The array must not be changed by you between a call of f01brf and a call of f04axf.

4:     $\mathbf{licn}$ – IntegerInput

On entry: the dimension of the arrays a and icn as declared in the (sub)program from which f01brf is called. Since the factorization is returned in a and icn, licn should be large enough to accommodate this and should ordinarily be $2$ to $4$ times as large as nz.

Constraint: $\mathbf{licn}\ge \mathbf{nz}$.

5:     $\mathbf{irn}\left(\mathbf{lirn}\right)$ – Integer arrayInput/Output

On entry: $\mathbf{irn}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{nz}$, must contain the row index of the nonzero element stored in $\mathbf{a}\left(i\right)$.

On exit: irn is overwritten and is not needed for subsequent calls of f01bsf or f04axf.

6:     $\mathbf{lirn}$ – IntegerInput

On entry: the dimension of the array irn as declared in the (sub)program from which f01brf is called. It need not be as large as licn; normally it will not need to be very much greater than nz.

Constraint: $\mathbf{lirn}\ge \mathbf{nz}$.

7:     $\mathbf{icn}\left(\mathbf{licn}\right)$ – Integer arrayCommunication Array

$\mathbf{icn}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{nz}$, must contain, on entry, the column index of the nonzero element stored in $\mathbf{a}\left(i\right)$. icn contains, on exit, the column indices of the nonzero elements in the factorization. The array must not be changed by you between a call of f01brf and subsequent calls of f01bsf or f04axf.

8:     $\mathbf{pivot}$ – Real (Kind=nag_wp)Input

On entry: should have a value in the range $0.0\le \mathbf{pivot}\le 0.9999$ and is used to control the choice of pivots. If $\mathbf{pivot}<0.0$, the value $0.0$ is assumed, and if $\mathbf{pivot}>0.9999$, the value $0.9999$ is assumed. When searching a row for a pivot, any element is excluded which is less than pivot times the largest of those elements in the row available as pivots. Thus decreasing pivot biases the algorithm to maintaining sparsity at the expense of stability.

Suggested value: $\mathbf{pivot}=0.1$ has been found to work well on test examples.

9:     $\mathbf{ikeep}\left(5\times \mathbf{n}\right)$ – Integer arrayCommunication Array

On exit: indexing information about the factorization.

You must not change ikeep between a call of f01brf and subsequent calls to f01bsf or f04axf.

10:   $\mathbf{iw}\left(8\times \mathbf{n}\right)$ – Integer arrayWorkspace

11:   $\mathbf{w}\left(\mathbf{n}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: if $\mathbf{grow}=\mathrm{.TRUE.}$, $\mathbf{w}\left(1\right)$ contains an estimate (an upper bound) of the increase in size of elements encountered during the factorization (see grow); the rest of the array is used as workspace.

If $\mathbf{grow}=\mathrm{.FALSE.}$, the array is not used.

12:   $\mathbf{lblock}$ – LogicalInput

On entry: if $\mathbf{lblock}=\mathrm{.TRUE.}$, the matrix is preordered into block lower triangular form before the $LU$ factorization is performed; otherwise the entire matrix is factorized.

Suggested value: $\mathbf{lblock}=\mathrm{.TRUE.}$ unless the matrix is known to be irreducible, or is singular and an upper bound on the rank is required.

13:   $\mathbf{grow}$ – LogicalInput

On entry: if $\mathbf{grow}=\mathrm{.TRUE.}$, then on exit $\mathbf{w}\left(1\right)$ contains an estimate (an upper bound) of the increase in size of elements encountered during the factorization. If the matrix is well-scaled (see Section 9.2), then a high value for $\mathbf{w}\left(1\right)$ indicates that the $LU$ factorization may be inaccurate and you should be wary of the results and perhaps increase the argument pivot for subsequent runs (see Section 7).

Suggested value: $\mathbf{grow}=\mathrm{.TRUE.}$.

14:   $\mathbf{abort}\left(4\right)$ – Logical arrayInput

On entry: if $\mathbf{abort}\left(1\right)=\mathrm{.TRUE.}$, f01brf will exit immediately on detecting a structural singularity (one that depends on the pattern of nonzeros) and return $\mathbf{ifail}=\mathbf{1}$; otherwise it will complete the factorization (see Section 9.3).

If $\mathbf{abort}\left(2\right)=\mathrm{.TRUE.}$, f01brf will exit immediately on detecting a numerical singularity (one that depends on the numerical values) and return $\mathbf{ifail}=\mathbf{2}$; otherwise it will complete the factorization (see Section 9.3).

If $\mathbf{abort}\left(3\right)=\mathrm{.TRUE.}$, f01brf will exit immediately (with $\mathbf{ifail}=\mathbf{5}$) when the arrays a and icn are filled up by the previously factorized, active and unfactorized parts of the matrix; otherwise it continues so that better guidance on necessary array sizes can be given in $\mathbf{idisp}\left(6\right)$ and $\mathbf{idisp}\left(7\right)$, and will exit with ifail in the range $4$ to $6$. Note that there is always an immediate error exit if the array irn is too small.

If $\mathbf{abort}\left(4\right)=\mathrm{.TRUE.}$, f01brf exits immediately (with $\mathbf{ifail}=\mathbf{13}$) if it finds duplicate elements in the input matrix.

If $\mathbf{abort}\left(4\right)=\mathrm{.FALSE.}$, f01brf proceeds using a value equal to the sum of the duplicate elements. In either case details of each duplicate element are output on the current advisory message unit (see x04abf), unless suppressed by the value of ifail on entry.

Suggested values:

• $\mathbf{abort}\left(1\right)=\mathrm{.TRUE.}$;
• $\mathbf{abort}\left(2\right)=\mathrm{.TRUE.}$;
• $\mathbf{abort}\left(3\right)=\mathrm{.FALSE.}$;
• $\mathbf{abort}\left(4\right)=\mathrm{.TRUE.}$.

15:   $\mathbf{idisp}\left(10\right)$ – Integer arrayCommunication Array

On exit: contains information about the factorization.

$\mathbf{idisp}\left(1\right)$ and $\mathbf{idisp}\left(2\right)$ indicate the position in arrays a and icn of the first and last elements in the $LU$ factorization of the diagonal blocks. ($\mathbf{idisp}\left(2\right)$ gives the number of nonzeros in the factorization.) $\mathbf{idisp}\left(1\right)$ and $\mathbf{idisp}\left(2\right)$ must not be changed by you between a call of f01brf and subsequent calls to f01bsf or f04axf.

$\mathbf{idisp}\left(3\right)$ and $\mathbf{idisp}\left(4\right)$ monitor the adequacy of ‘elbow room’ in the arrays irn and a (and icn) respectively, by giving the number of times that the data in these arrays has been compressed during the factorization to release more storage. If either $\mathbf{idisp}\left(3\right)$ or $\mathbf{idisp}\left(4\right)$ is quite large (say greater than $10$), it will probably pay you to increase the size of the corresponding array(s) for subsequent runs. If either is very low or zero, then you can perhaps save storage by reducing the size of the corresponding array(s).

$\mathbf{idisp}\left(5\right)$, when $\mathbf{lblock}=\mathrm{.FALSE.}$, gives an upper bound on the rank of the matrix; when $\mathbf{lblock}=\mathrm{.TRUE.}$, gives an upper bound on the sum of the ranks of the lower triangular blocks.

$\mathbf{idisp}\left(6\right)$ and $\mathbf{idisp}\left(7\right)$ give the minimum size of arrays irn and a (and icn) respectively which would enable a successful run on an identical matrix (but some ‘elbow-room’ should be allowed – see Section 9).

$\mathbf{idisp}\left(8\right)$ to $\left(10\right)$ are only used if $\mathbf{lblock}=\mathrm{.TRUE.}$.

• $\mathbf{idisp}\left(8\right)$ gives the structural rank of the matrix.
• $\mathbf{idisp}\left(9\right)$ gives the number of diagonal blocks.
• $\mathbf{idisp}\left(10\right)$ gives the size of the largest diagonal block.

You must not change idisp between a call of f01brf and subsequent calls to f01bsf or f04axf.

16:   $\mathbf{ifail}$ – IntegerInput/Output

For this routine, the normal use of ifail is extended to control the printing of error and warning messages as well as specifying hard or soft failure (see Section 3.4 in How to Use the NAG Library and its Documentation).

On entry: ifail must be set to a value with the decimal expansion $\mathit{cba}$, where each of the decimal digits $c$, $b$ and $a$ must have a value of $0$ or $1$.

$a=0$	specifies hard failure, otherwise soft failure;
$b=0$	suppresses error messages, otherwise error messages will be printed (see Section 6);
$c=0$	suppresses warning messages, otherwise warning messages will be printed (see Section 6).

The recommended value for inexperienced users is $110$ (i.e., hard failure with all messages printed).

On exit: $\mathbf{ifail}=\mathbf{0}$ unless the routine detects an error or a warning has been flagged (see Section 6).

6

Error Indicators and Warnings

$\mathbf{ifail}=-2$

Successful factorization of a numerically singular matrix (which may also be structurally singular) (see Section 9.3).

$\mathbf{ifail}=-1$

Successful factorization of a structurally singular matrix (see Section 9.3).

$\mathbf{ifail}=1$

The matrix is structurally singular and the factorization has been abandoned ($\mathbf{abort}\left(1\right)$ was .TRUE. on entry).

$\mathbf{ifail}=2$

The matrix is numerically singular and the factorization has been abandoned ($\mathbf{abort}\left(2\right)$ was .TRUE. on entry).

$\mathbf{ifail}=3$

lirn is too small: there is not enough space in the array irn to continue the factorization. You are recommended to try again with lirn (and the length of irn) equal to at least $\mathbf{idisp}\left(6\right)+\mathbf{n}/2$.

$\mathbf{ifail}=4$

licn is much too small: there is much too little space in the arrays a and icn to continue the factorization.

$\mathbf{ifail}=5$

licn is too small: there is not enough space in the arrays a and icn to store the factorization. If $\mathbf{abort}\left(3\right)$ was .FALSE. on entry, the factorization has been completed but some of the $LU$ factors have been discarded to create space; $\mathbf{idisp}\left(7\right)$ then gives the minimum value of licn (i.e., the minimum length of a and icn) required for a successful factorization of the same matrix.

$\mathbf{ifail}=6$

licn and lirn are both too small: effectively this is a combination of $\mathbf{ifail}=\mathbf{3}$ and $\mathbf{5}$ (with $\mathbf{abort}\left(3\right)=\mathrm{.FALSE.}$).

$\mathbf{ifail}=7$

licn is too small: there is not enough space in the arrays a and icn for the permutation to block triangular form.

$\mathbf{ifail}=8$

On entry,

$\mathbf{n}\le 0$.

$\mathbf{ifail}=9$

On entry,

$\mathbf{nz}\le 0$.

$\mathbf{ifail}=10$

On entry,

$\mathbf{licn}<\mathbf{nz}$.

$\mathbf{ifail}=11$

On entry,

$\mathbf{lirn}<\mathbf{nz}$.

$\mathbf{ifail}=12$

On entry, an element of the input matrix has a row or column index (i.e., an element of irn or icn) outside the range $1$ to n.

$\mathbf{ifail}=13$

Duplicate elements have been found in the input matrix and the factorization has been abandoned ($\mathbf{abort}\left(4\right)=\mathrm{.TRUE.}$ on entry).

$\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$\mathbf{ifail}=-999$

Dynamic memory allocation failed.

See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7

Accuracy

8

Parallelism and Performance

9

Further Comments

9.1

Timing

9.2

Scaling

9.3

Singular and Rectangular Systems

9.4

Duplicated Nonzeros

9.5

Determinant

   det = 1.0
   id = idisp(1)
   Do 10 i = 1, n
      idg = id + ikeep(3*n+i)
      det = det*a(idg)
      If (ikeep(n+i).ne.i)det = -det
      If (ikeep(2*n+i).ne.i)det = -det
      id = id + ikeep(i)
10 Continue

10

Example

10.1

Program Text

Program Text (f01brfe.f90)

10.2

Program Data

Program Data (f01brfe.d)

10.3

Program Results

Program Results (f01brfe.r)