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$ 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$, nag_sparse_direct_real_gen_lu (f11me) computes an $LU$ factorization of $A$ with partial pivoting, ${P}_{r}A{P}_{c}=LU$, where ${P}_{r}$ is a row permutation matrix (computed by nag_sparse_direct_real_gen_lu (f11me)), ${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 nag_sparse_direct_real_gen_setup (f11md). 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 arguments (order and number of nonzero values) of the matrix $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:     $\mathrm{n}$int64int32nag_int scalar
$n$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
2:     $\mathrm{irowix}\left(:\right)$int64int32nag_int array
The dimension of the array irowix must be at least $\mathit{nnz}$, the number of nonzeros of the sparse matrix $A$
The row index array of sparse matrix $A$.
3:     $\mathrm{a}\left(:\right)$ – double array
The dimension of the array a must be at least $\mathit{nnz}$, the number of nonzeros of the sparse matrix $A$
The array of nonzero values in the sparse matrix $A$.
4:     $\mathrm{iprm}\left(7×{\mathbf{n}}\right)$int64int32nag_int array
Contains the column permutation which defines the permutation ${P}_{c}$ and associated data structures as computed by function nag_sparse_direct_real_gen_setup (f11md).
5:     $\mathrm{thresh}$ – double scalar
Suggested value: ${\mathbf{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 ${\mathbf{thresh}}=0.0$ if you are confident about the stability of the factorization, for example, if $A$ is diagonally dominant.
The diagonal pivoting threshold, $t$. At step $j$ of the Gaussian elimination, if $\left|{A}_{jj}\right|\ge t\left(\underset{i\ge j}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left|{A}_{ij}\right|\right)$, use ${A}_{jj}$ as a pivot, otherwise use $\underset{i\ge j}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\left|{A}_{ij}\right|$. 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).
Constraint: $0.0\le {\mathbf{thresh}}\le 1.0$.
6:     $\mathrm{nzlmx}$int64int32nag_int scalar
Indicates the available size of array il. The dimension of il should be at least $7×{\mathbf{n}}+{\mathbf{nzlmx}}+4$. A good range for nzlmx that works for many problems is $\mathit{nnz}$ to $8×\mathit{nnz}$, where $\mathit{nnz}$ is the number of nonzeros in the sparse matrix $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: ${\mathbf{nzlmx}}\ge 1$.
7:     $\mathrm{nzlumx}$int64int32nag_int scalar
Indicates the available size of array lval. The dimension of lval should be at least nzlumx.
Constraint: ${\mathbf{nzlumx}}\ge 1$.
8:     $\mathrm{nzumx}$int64int32nag_int scalar
Indicates the available sizes of arrays iu and uval. The dimension of iu should be at least $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 $\mathit{nnz}$ to $8×\mathit{nnz}$, where $\mathit{nnz}$ is the number of nonzeros in the sparse matrix $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: ${\mathbf{nzumx}}\ge 1$.

None.

### Output Parameters

1:     $\mathrm{iprm}\left(7×{\mathbf{n}}\right)$int64int32nag_int array
Part of the array is modified to record the row permutation ${P}_{r}$ determined by pivoting.
2:     $\mathrm{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:     $\mathrm{il}\left(7×{\mathbf{n}}+{\mathbf{nzlmx}}+4\right)$int64int32nag_int array
Encapsulates the sparsity pattern of matrix $L$.
4:     $\mathrm{lval}\left({\mathbf{nzlumx}}\right)$ – double array
Records the nonzero values of matrix $L$ and some of the nonzero values of matrix $U$.
5:     $\mathrm{iu}\left(2×{\mathbf{n}}+{\mathbf{nzumx}}+1\right)$int64int32nag_int array
Encapsulates the sparsity pattern of matrix $U$.
6:     $\mathrm{uval}\left({\mathbf{nzumx}}\right)$ – double array
Records some of the nonzero values of matrix $U$.
7:     $\mathrm{nnzl}$int64int32nag_int scalar
The number of nonzero values in the matrix $L$.
8:     $\mathrm{nnzu}$int64int32nag_int scalar
The number of nonzero values in the matrix $U$.
9:     $\mathrm{flop}$ – double scalar
The number of floating-point operations performed.
10:   $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\mathbf{ifail}}=1$
Constraint: $0.0\le {\mathbf{thresh}}\le 1.0$.
Constraint: ${\mathbf{n}}\ge 0$.
Constraint: ${\mathbf{nzlmx}}\ge 1$.
Constraint: ${\mathbf{nzlumx}}\ge 1$.
Constraint: ${\mathbf{nzumx}}\ge 1$.
${\mathbf{ifail}}=2$
Insufficient nzlmx.
${\mathbf{ifail}}=3$
Insufficient nzumx.
${\mathbf{ifail}}=4$
Insufficient nzlumx.
${\mathbf{ifail}}=5$
The matrix is singular – no factorization possible.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

The computed factors $L$ and $U$ are the exact factors of a perturbed matrix $A+E$, where
 $E ≤ c n ε L U ,$
$c\left(n\right)$ is a modest linear function of $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 call to nag_sparse_direct_real_gen_lu (f11me) may be followed by calls to the functions:

## Example

This example computes the $LU$ factorization of the matrix $A$, where
 $A= 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 .$
```function f11me_example

fprintf('f11me example results\n\n');

% LU factorize sparse A and solve Ax = b
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];

% Permute matrix
iprm   = zeros(1, 7*n, 'int64');
spec   = 'M';
% Calculate COLAMD permutation
[iprm, ifail] = f11md( ...
spec, n, icolzp, irowix, iprm);

% Factorize
thresh = 1;
nzlmx  = int64(8*nz);
nzlumx = nzlmx;
nzumx  = nzlmx;
[iprm, nzlumx, il, lval, iu, uval, nnzl, nnzu, flop, ifail] = ...
f11me(...
n, irowix, a, iprm, thresh, nzlmx, nzlumx, nzumx);

fprintf('\nNumber of nonzeros in factors (excluding unit diagonal)');
fprintf('%8d\n',nnzl + nnzu - n);

disp('Factor elements in lval');
fprintf('%7.2f%7.2f%7.2f%7.2f%7.2f\n',lval(1:10)');
disp('Factor elements in uval');
fprintf('%7.2f%7.2f%7.2f%7.2f\n',uval(1:4)');

```
```f11me example results

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