PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_matop_real_gen_sparse_lu_reuse (f01bs)
Purpose
nag_matop_real_gen_sparse_lu_reuse (f01bs) factorizes a real sparse matrix using the pivotal sequence previously obtained by
nag_matop_real_gen_sparse_lu (f01br) when a matrix of the same sparsity pattern was factorized.
Syntax
[
a,
w,
rpmin,
ifail] = f01bs(
n,
a,
ivect,
jvect,
icn,
ikeep,
grow,
idisp, 'nz',
nz, 'licn',
licn, 'eta',
eta, 'abort',
abort)
[
a,
w,
rpmin,
ifail] = nag_matop_real_gen_sparse_lu_reuse(
n,
a,
ivect,
jvect,
icn,
ikeep,
grow,
idisp, 'nz',
nz, 'licn',
licn, 'eta',
eta, 'abort',
abort)
Description
nag_matop_real_gen_sparse_lu_reuse (f01bs) accepts as input a real sparse matrix of the same sparsity pattern as a matrix previously factorized by a call of
nag_matop_real_gen_sparse_lu (f01br). It first applies to the matrix the same permutations as were used by
nag_matop_real_gen_sparse_lu (f01br), both for permutation to block triangular form and for pivoting, and then performs Gaussian elimination to obtain the
$LU$ factorization of the diagonal blocks.
Extensive data checks are made; duplicated nonzeros can be accumulated.
The factorization is intended to be used by
nag_linsys_real_sparse_fac_solve (f04ax) to solve sparse systems of linear equations
$Ax=b$ or
${A}^{\mathrm{T}}x=b$.
nag_matop_real_gen_sparse_lu_reuse (f01bs) is much faster than
nag_matop_real_gen_sparse_lu (f01br) and in some applications it is expected that there will be many calls of
nag_matop_real_gen_sparse_lu_reuse (f01bs) for each call of
nag_matop_real_gen_sparse_lu (f01br).
The method is fully described in
Duff (1977).
A more recent algorithm for the same calculation is provided by
nag_sparse_direct_real_gen_lu (f11me).
References
Duff I S (1977) MA28 – a set of Fortran subroutines for sparse unsymmetric linear equations AERE Report R8730 HMSO
Parameters
Compulsory Input Parameters
 1:
$\mathrm{n}$ – int64int32nag_int scalar

$n$, the order of the matrix $A$.
Constraint:
${\mathbf{n}}>0$.
 2:
$\mathrm{a}\left({\mathbf{licn}}\right)$ – double array

${\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 nag_matop_real_gen_sparse_lu_reuse (f01bs) will reorder them.
 3:
$\mathrm{ivect}\left({\mathbf{nz}}\right)$ – int64int32nag_int array
 4:
$\mathrm{jvect}\left({\mathbf{nz}}\right)$ – int64int32nag_int array

${\mathbf{ivect}}\left(\mathit{i}\right)$ and ${\mathbf{jvect}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{nz}}$, must contain the row index and the column index respectively of the nonzero element stored in ${\mathbf{a}}\left(i\right)$.
 5:
$\mathrm{icn}\left({\mathbf{licn}}\right)$ – int64int32nag_int array

icn contains, on entry, the same information as output by
nag_matop_real_gen_sparse_lu (f01br). It must not be changed by you between a call of
nag_matop_real_gen_sparse_lu_reuse (f01bs) and a call of
nag_linsys_real_sparse_fac_solve (f04ax).
icn is used as internal workspace prior to being restored on exit and hence is unchanged.
 6:
$\mathrm{ikeep}\left(5\times {\mathbf{n}}\right)$ – int64int32nag_int array

The same indexing information about the factorization as output in
ikeep by
nag_matop_real_gen_sparse_lu (f01br).
You must
not change
ikeep between a call of
nag_matop_real_gen_sparse_lu_reuse (f01bs) and subsequent calls to
nag_linsys_real_sparse_fac_solve (f04ax).
 7:
$\mathrm{grow}$ – logical scalar

If
${\mathbf{grow}}=\mathit{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 wellscaled (see
Further Comments), 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
Accuracy).
 8:
$\mathrm{idisp}\left(2\right)$ – int64int32nag_int array

${\mathbf{idisp}}\left(1\right)$ and
${\mathbf{idisp}}\left(2\right)$ must be as output in
idisp by the previous call of
nag_matop_real_gen_sparse_lu (f01br).
Optional Input Parameters
 1:
$\mathrm{nz}$ – int64int32nag_int scalar

Default:
the dimension of the arrays
ivect,
jvect. (An error is raised if these dimensions are not equal.)
The number of nonzero elements in the matrix $A$.
Constraint:
${\mathbf{nz}}>0$.
 2:
$\mathrm{licn}$ – int64int32nag_int scalar

Default:
the dimension of the arrays
a,
icn. (An error is raised if these dimensions are not equal.)
The dimension of the arrays
a and
icn. it should have the same value as it had for
nag_matop_real_gen_sparse_lu (f01br).
Constraint:
${\mathbf{licn}}\ge {\mathbf{nz}}$.
 3:
$\mathrm{eta}$ – double scalar
Default:
$0.0001$
The relative pivot threshold below which an error diagnostic is provoked and
ifail is set to
${\mathbf{ifail}}={\mathbf{7}}$. If
eta is greater than
$1.0$, then no check on pivot size is made.
 4:
$\mathrm{abort}$ – logical scalar
Default:
$\mathit{true}$
If
${\mathbf{abort}}=\mathit{true}$,
nag_matop_real_gen_sparse_lu_reuse (f01bs) exits immediately (with
${\mathbf{ifail}}={\mathbf{8}}$) if it finds duplicate elements in the input matrix.
If ${\mathbf{abort}}=\mathit{false}$, nag_matop_real_gen_sparse_lu_reuse (f01bs) 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
nag_file_set_unit_advisory (x04ab)), unless suppressed by the value of
ifail on entry.
Output Parameters
 1:
$\mathrm{a}\left({\mathbf{licn}}\right)$ – double array

The nonzero elements in the
$LU$ factorization. The array must
not be changed by you between a call of
nag_matop_real_gen_sparse_lu_reuse (f01bs) and a call of
nag_linsys_real_sparse_fac_solve (f04ax).
 2:
$\mathrm{w}\left({\mathbf{n}}\right)$ – double array

If
${\mathbf{grow}}=\mathit{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}}=\mathit{false}$, the array is not used.
 3:
$\mathrm{rpmin}$ – double scalar

If
eta is less than
$1.0$, then
rpmin gives the smallest ratio of the pivot to the largest element in the row of the corresponding upper triangular factor thus monitoring the stability of the factorization. If
rpmin is very small it may be advisable to perform a new factorization using
nag_matop_real_gen_sparse_lu (f01br).
 4:
$\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:
Cases prefixed with W are classified as warnings and
do not generate an error of type NAG:error_n. See nag_issue_warnings.
 ${\mathbf{ifail}}=1$

On entry,  ${\mathbf{n}}\le 0$. 
 ${\mathbf{ifail}}=2$

On entry,  ${\mathbf{nz}}\le 0$. 
 ${\mathbf{ifail}}=3$

On entry,  ${\mathbf{licn}}<{\mathbf{nz}}$. 
 ${\mathbf{ifail}}=4$

On entry, an element of the input matrix has a row or column index (i.e., an element of
ivect or
jvect) outside the range
$1$ to
n.
 ${\mathbf{ifail}}=5$

The input matrix is incompatible with the matrix factorized by the previous call of
nag_matop_real_gen_sparse_lu (f01br) (see
Further Comments).
 ${\mathbf{ifail}}=6$

The input matrix is numerically singular.
 W ${\mathbf{ifail}}=7$

A very small pivot has been detected (see
Arguments,
eta). The factorization has been completed but is potentially unstable.
 ${\mathbf{ifail}}=8$

Duplicate elements have been found in the input matrix and the factorization has been abandoned (${\mathbf{abort}}=\mathit{true}$ on entry).
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
 ${\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 factorization obtained is exact for a perturbed matrix whose
$\left(i,j\right)$th element differs from
${a}_{ij}$ by less than
$3\epsilon \rho {m}_{ij}$ where
$\epsilon $ is the
machine precision,
$\rho $ is the growth value returned in
${\mathbf{w}}\left(1\right)$ if
${\mathbf{grow}}=\mathit{true}$, and
${m}_{ij}$ the number of Gaussian elimination operations applied to element
$\left(i,j\right)$.
If
$\rho ={\mathbf{w}}\left(1\right)$ is very large or
rpmin is very small, then a fresh call of
nag_matop_real_gen_sparse_lu (f01br) is recommended.
Further Comments
If you have a sequence of problems with the same sparsity pattern then
nag_matop_real_gen_sparse_lu_reuse (f01bs) is recommended after
nag_matop_real_gen_sparse_lu (f01br) has been called for one such problem. It is typically
$4$ to
$7$ times faster but is potentially unstable since the previous pivotal sequence is used. Further details on timing are given in the document for
nag_matop_real_gen_sparse_lu (f01br).
If growth estimation is performed (${\mathbf{grow}}=\mathit{true}$), then the time increases by between $5\%$ and $10\%$. Pivot size monitoring (${\mathbf{eta}}\le 1.0$) involves a similar overhead.
We normally expect this function to be entered with a matrix having the same pattern of nonzeros as was earlier presented to
nag_matop_real_gen_sparse_lu (f01br). However there is no record of this pattern, but rather a record of the pattern including all fillins. Therefore we permit additional nonzeros in positions corresponding to fillins.
If singular matrices are being treated then it is also required that the present matrix be sufficiently like the previous one for the same permutations to be suitable for factorization with the same set of zero pivots.
Example
This example factorizes the real sparse matrices
and
This example program simply prints the values of
${\mathbf{w}}\left(1\right)$ and
rpmin returned by
nag_matop_real_gen_sparse_lu_reuse (f01bs). Normally the calls of
nag_matop_real_gen_sparse_lu (f01br) and
nag_matop_real_gen_sparse_lu_reuse (f01bs) would be followed by calls of
nag_linsys_real_sparse_fac_solve (f04ax).
Open in the MATLAB editor:
f01bs_example
function f01bs_example
fprintf('f01bs example results\n\n');
n = int64(6);
nz = int64(15);
z = int64(15);
a = zeros(150,1);
a(1:15) = [5; 2;1;2; 3; 2;1;1; 1;1;2;3; 1;1;6];
irn = zeros(75,1,'int64');
icn = zeros(150,1,'int64');
irn(1:15) = [int64(1); 2;2;2; 3; 4;4;4; 5;5;5;5; 6;6;6];
icn(1:15) = [int64(1); 2;3;4; 3; 1;4;5; 1;4;5;6; 1;2;6];
abort = [true; true; false; true];
ivect = irn(1:15);
jvect = icn(1:15);
grow = true;
[a, irn, icn, ikeep, w, idisp, ifail] = ...
f01br( ...
n, nz, a, irn, icn, abort);
fprintf('For first matrix:\nValue of w(1) = %7.4f\n\n', w(1));
a(1:15) = [10; 12;3;1; 15; 2;10;1; 1;5;1;1; 1;2;6 ];
eta = 0.1;
[a, w, rpmin, ifail] = f01bs( ...
n, a, ivect, jvect, icn, ikeep, grow, idisp, ...
'eta',eta);
fprintf('For second matrix:\nValue of w(1) = %7.4f\n\n', w(1));
fprintf('Value of rpmin = %7.4f\n', rpmin);
f01bs example results
For first matrix:
Value of w(1) = 18.0000
For second matrix:
Value of w(1) = 51.0000
Value of rpmin = 0.1000
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015