Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_sparse_direct_real_gen_lu (f11me)

## Purpose

nag_sparse_direct_real_gen_lu (f11me) computes the LU $LU$ factorization of a real sparse matrix in compressed column (Harwell–Boeing), column-permuted format.

## Syntax

[iprm, nzlumx, il, lval, iu, uval, nnzl, nnzu, flop, ifail] = f11me(n, irowix, a, iprm, thresh, nzlmx, nzlumx, nzumx)
[iprm, nzlumx, il, lval, iu, uval, nnzl, nnzu, flop, ifail] = nag_sparse_direct_real_gen_lu(n, irowix, a, iprm, thresh, nzlmx, nzlumx, nzumx)

## Description

Given a real sparse matrix A$A$, nag_sparse_direct_real_gen_lu (f11me) computes an LU $LU$ factorization of A$A$ with partial pivoting, Pr A Pc = LU ${P}_{r}A{P}_{c}=LU$, where Pr${P}_{r}$ is a row permutation matrix (computed by nag_sparse_direct_real_gen_lu (f11me)), Pc${P}_{c}$ is a (supplied) column permutation matrix, L$L$ is unit lower triangular and U$U$ is upper triangular. The column permutation matrix, Pc${P}_{c}$, must be computed by a prior call to nag_sparse_direct_real_gen_setup (f11md). The matrix A$A$ must be presented in the column permuted, compressed column (Harwell–Boeing) format.
The LU $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$L$ and U$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$A$. This is due to unpredictable fill-in created by partial pivoting. nag_sparse_direct_real_gen_lu (f11me) 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.

## 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

## Parameters

### Compulsory Input Parameters

1:     n – int64int32nag_int scalar
n$n$, the order of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.
2:     irowix( : $:$) – int64int32nag_int array
Note: the dimension of the array irowix must be at least nnz$\mathit{nnz}$, the number of nonzeros of the sparse matrix A$A$.
The row index array of sparse matrix A$A$.
3:     a( : $:$) – double array
Note: the dimension of the array a must be at least nnz$\mathit{nnz}$, the number of nonzeros of the sparse matrix A$A$.
The array of nonzero values in the sparse matrix A$A$.
4:     iprm(7 × n$7×{\mathbf{n}}$) – int64int32nag_int array
Contains the column permutation which defines the permutation Pc${P}_{c}$ and associated data structures as computed by function nag_sparse_direct_real_gen_setup (f11md).
5:     thresh – double scalar
The diagonal pivoting threshold, t$t$. At step j$j$ of the Gaussian elimination, if |Ajj|t(maxij |Aij|)$|{A}_{jj}|\ge t\left(\underset{i\ge j}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}|{A}_{ij}|\right)$, use Ajj${A}_{jj}$ as a pivot, otherwise use maxij |Aij|$\underset{i\ge j}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}|{A}_{ij}|$. A value of t = 1$t=1$ corresponds to partial pivoting, a value of t = 0$t=0$ corresponds to always choosing the pivot on the diagonal (unless it is zero).
Constraint: 0.0thresh1.0$0.0\le {\mathbf{thresh}}\le 1.0$.
6:     nzlmx – int64int32nag_int scalar
Indicates the available size of array il. The dimension of il should be at least 7 × n + nzlmx + 4$7×{\mathbf{n}}+{\mathbf{nzlmx}}+4$. A good range for nzlmx that works for many problems is nnz$\mathit{nnz}$ to 8 × nnz$8×\mathit{nnz}$, where nnz$\mathit{nnz}$ is the number of nonzeros in the sparse matrix A$A$. If, on exit, ${\mathbf{ifail}}={\mathbf{2}}$, the given nzlmx was too small and you should attempt to provide more storage and call the function again.
Constraint: nzlmx1${\mathbf{nzlmx}}\ge 1$.
7:     nzlumx – int64int32nag_int scalar
Indicates the available size of array lval. The dimension of lval should be at least nzlumx.
Constraint: nzlumx1${\mathbf{nzlumx}}\ge 1$.
8:     nzumx – int64int32nag_int scalar
Indicates the available sizes of arrays iu and uval. The dimension of iu should be at least 2 × n + nzumx + 1$2×{\mathbf{n}}+{\mathbf{nzumx}}+1$ and the dimension of uval should be at least nzumx. A good range for nzumx that works for many problems is nnz$\mathit{nnz}$ to 8 × nnz$8×\mathit{nnz}$, where nnz$\mathit{nnz}$ is the number of nonzeros in the sparse matrix A$A$. If, on exit, ${\mathbf{ifail}}={\mathbf{3}}$, the given nzumx was too small and you should attempt to provide more storage and call the function again.
Constraint: nzumx1${\mathbf{nzumx}}\ge 1$.

None.

None.

### Output Parameters

1:     iprm(7 × n$7×{\mathbf{n}}$) – int64int32nag_int array
Part of the array is modified to record the row permutation Pr${P}_{r}$ determined by pivoting.
2:     nzlumx – int64int32nag_int scalar
If ${\mathbf{ifail}}={\mathbf{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 function again.
3:     il(7 × n + nzlmx + 4$7×{\mathbf{n}}+{\mathbf{nzlmx}}+4$ ) – int64int32nag_int array
Encapsulates the sparsity pattern of matrix L$L$.
4:     lval(nzlumx) – double array
Records the nonzero values of matrix L$L$ and some of the nonzero values of matrix U$U$.
5:     iu(2 × n + nzumx + 1$2×{\mathbf{n}}+{\mathbf{nzumx}}+1$) – int64int32nag_int array
Encapsulates the sparsity pattern of matrix U$U$.
6:     uval(nzumx${\mathbf{nzumx}}$) – double array
Records some of the nonzero values of matrix U$U$.
7:     nnzl – int64int32nag_int scalar
The number of nonzero values in the matrix L$L$.
8:     nnzu – int64int32nag_int scalar
The number of nonzero values in the matrix U$U$.
9:     flop – double scalar
The number of floating point operations performed.
10:   ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, n < 0${\mathbf{n}}<0$, or nzlmx < 1${\mathbf{nzlmx}}<1$, or nzlumx < 1${\mathbf{nzlumx}}<1$, or nzumx < 1${\mathbf{nzumx}}<1$, or thresh < 0.0${\mathbf{thresh}}<0.0$, or thresh > 1.0${\mathbf{thresh}}>1.0$.
ifail = 2${\mathbf{ifail}}=2$
nzlmx${\mathbf{nzlmx}}$ was not large enough. You should repeat the call with a larger value of nzlmx${\mathbf{nzlmx}}$, providing more storage for the output array il${\mathbf{il}}$.
ifail = 3${\mathbf{ifail}}=3$
nzumx${\mathbf{nzumx}}$ was not large enough. You should repeat the call with a larger value of nzumx${\mathbf{nzumx}}$, providing more storage for the output arrays iu${\mathbf{iu}}$ and uval${\mathbf{uval}}$.
ifail = 4${\mathbf{ifail}}=4$
nzlumx${\mathbf{nzlumx}}$ was not large enough. You should repeat the call with the value of nzlumx${\mathbf{nzlumx}}$ returned on exit, providing more storage for the output array lval${\mathbf{lval}}$.
ifail = 5${\mathbf{ifail}}=5$
The matrix A$A$ is singular and no factorization will be attempted.
ifail = 301${\mathbf{ifail}}=301$
Unable to allocate required internal workspace.

## Accuracy

The computed factors L$L$ and U$U$ are the exact factors of a perturbed matrix A + E$A+E$, where
 |E| ≤ c (n) ε |L| |U| , $| E | ≤ c ( n ) ε |L| |U| ,$
c(n)$c\left(n\right)$ is a modest linear function of n$n$, and ε$\epsilon$ 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 nag_sparse_direct_real_gen_diag (f11mm) after nag_sparse_direct_real_gen_lu (f11me) can help determine the quality of the factorization.

The total number of floating point operations depends on the sparsity pattern of the matrix A$A$.
A call to nag_sparse_direct_real_gen_lu (f11me) may be followed by calls to the functions:

## Example

```function nag_sparse_direct_real_gen_lu_example
n = int64(5);
nz = int64(11);
icolzp = [int64(1); 3; 5; 7; 9; 12];
irowix = [int64(1); 3; 1; 5; 2; 3; 2; 4; 3; 4; 5];
a = [2; 4; 1; -2; 1; 1; -1; 1; 1; 2; 3];
iprm = zeros(1, 7*n, 'int64');
spec = 'M';
thresh = 1;
nzlmx = int64(8*nz);
nzlumx = int64(8*nz);
nzumx = int64(8*nz);

% Calculate COLAMD permutation
[iprm, ifail] = ...
nag_sparse_direct_real_gen_setup(spec, n, icolzp, irowix, iprm);

% Factorise
[iprm, nzlumx, il, lval, iu, uval, nnzl, nnzu, flop, ifail] = ...
nag_sparse_direct_real_gen_lu(n, irowix, a, iprm, thresh, nzlmx, nzlumx, nzumx);
iprm, nzlumx, nnzl, nnzu, flop, ifail

fprintf('\nNumber of nonzeros in factors (excluding unit diagonal)\n%10d\n', nnzl+nnzu-n);
fprintf('\nFactor elements in lval\n');
fprintf('%6.2f',lval(1:10));
fprintf('\n');
fprintf('\nFactor elements in uval\n');
fprintf('%6.2f',uval(1:4));
fprintf('\n');
```
```

iprm =

Columns 1 through 4

1                    0                    4                    3

Columns 5 through 8

2                    4                    3                    1

Columns 9 through 12

2                    0                    2                    0

Columns 13 through 16

8                    6                    4                    4

Columns 17 through 20

2                   11                    8                    6

Columns 21 through 24

1                    2                    3                    4

Columns 25 through 28

5                    2                    2                    2

Columns 29 through 32

2                    1                    1                    1

Columns 33 through 35

1                    2                    0

nzlumx =

88

nnzl =

9

nnzu =

10

flop =

19

ifail =

0

Number of nonzeros in factors (excluding unit diagonal)
14

Factor elements in lval
-2.00 -0.50  4.00  0.50  2.00  0.50 -1.00  0.50  1.00 -1.00

Factor elements in uval
1.00  3.00  1.00  1.00

```
```function f11me_example
n = int64(5);
nz = int64(11);
icolzp = [int64(1); 3; 5; 7; 9; 12];
irowix = [int64(1); 3; 1; 5; 2; 3; 2; 4; 3; 4; 5];
a = [2; 4; 1; -2; 1; 1; -1; 1; 1; 2; 3];
iprm = zeros(1, 7*n, 'int64');
spec = 'M';
thresh = 1;
nzlmx = int64(8*nz);
nzlumx = int64(8*nz);
nzumx = int64(8*nz);

% Calculate COLAMD permutation
[iprm, ifail] = f11md(spec, n, icolzp, irowix, iprm);

% Factorise
[iprm, nzlumx, il, lval, iu, uval, nnzl, nnzu, flop, ifail] = ...
f11me(n, irowix, a, iprm, thresh, nzlmx, nzlumx, nzumx);
iprm, nzlumx, nnzl, nnzu, flop, ifail

fprintf('\nNumber of nonzeros in factors (excluding unit diagonal)\n%10d\n', nnzl+nnzu-n);
fprintf('\nFactor elements in lval\n');
fprintf('%6.2f',lval(1:10));
fprintf('\n');
fprintf('\nFactor elements in uval\n');
fprintf('%6.2f',uval(1:4));
fprintf('\n');
```
```

iprm =

Columns 1 through 4

1                    0                    4                    3

Columns 5 through 8

2                    4                    3                    1

Columns 9 through 12

2                    0                    2                    0

Columns 13 through 16

8                    6                    4                    4

Columns 17 through 20

2                   11                    8                    6

Columns 21 through 24

1                    2                    3                    4

Columns 25 through 28

5                    2                    2                    2

Columns 29 through 32

2                    1                    1                    1

Columns 33 through 35

1                    2                    0

nzlumx =

88

nnzl =

9

nnzu =

10

flop =

19

ifail =

0

Number of nonzeros in factors (excluding unit diagonal)
14

Factor elements in lval
-2.00 -0.50  4.00  0.50  2.00  0.50 -1.00  0.50  1.00 -1.00

Factor elements in uval
1.00  3.00  1.00  1.00

```