hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

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 LULU 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 = bAx=b or ATx = bATx=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:     n – int64int32nag_int scalar
nn, the order of the matrix AA.
Constraint: n > 0n>0.
2:     a(licn) – double array
licn, the dimension of the array, must satisfy the constraint licnnzlicnnz.
a(i)ai, for i = 1,2,,nzi=1,2,,nz, must contain the nonzero elements of the sparse matrix AA. They can be in any order since nag_matop_real_gen_sparse_lu_reuse (f01bs) will reorder them.
3:     ivect(nz) – int64int32nag_int array
4:     jvect(nz) – int64int32nag_int array
nz, the dimension of the array, must satisfy the constraint nz > 0nz>0.
ivect(i)ivecti and jvect(i)jvecti, for i = 1,2,,nzi=1,2,,nz, must contain the row index and the column index respectively of the nonzero element stored in a(i)ai.
5:     icn(licn) – int64int32nag_int array
licn, the dimension of the array, must satisfy the constraint licnnzlicnnz.
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:     ikeep(5 × n5×n) – 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:     grow – logical scalar
If grow = truegrow=true, then on exit w(1)w1 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 [Further Comments]), then a high value for w(1)w1 indicates that the LULU factorization may be inaccurate and you should be wary of the results and perhaps increase the parameter pivot for subsequent runs (see Section [Accuracy]).
8:     idisp(22) – int64int32nag_int array
idisp(1)idisp1 and idisp(2)idisp2 must be as output in idisp by the previous call of nag_matop_real_gen_sparse_lu (f01br).

Optional Input Parameters

1:     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 AA.
Constraint: nz > 0nz>0.
2:     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 as declared in the (sub)program from which nag_matop_real_gen_sparse_lu_reuse (f01bs) is called. It should have the same value as it had for nag_matop_real_gen_sparse_lu (f01br).
Constraint: licnnzlicnnz.
3:     eta – double scalar
The relative pivot threshold below which an error diagnostic is provoked and ifail is set to ifail = 7ifail=7. If eta is greater than 1.01.0, then no check on pivot size is made.
Default: 0.00010.0001
4:     abort – logical scalar
If abort = trueabort=true, nag_matop_real_gen_sparse_lu_reuse (f01bs) exits immediately (with ifail = 8ifail=8) if it finds duplicate elements in the input matrix.
If abort = falseabort=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.
Default: truetrue

Input Parameters Omitted from the MATLAB Interface

iw

Output Parameters

1:     a(licn) – double array
The nonzero elements in the LULU 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:     w(n) – double array
If grow = truegrow=true, w(1)w1 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 grow = falsegrow=false, the array is not used.
3:     rpmin – double scalar
If eta is less than 1.01.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:     ifail – int64int32nag_int scalar
ifail = 0ifail=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.

  ifail = 1ifail=1
On entry,n0n0.
  ifail = 2ifail=2
On entry,nz0nz0.
  ifail = 3ifail=3
On entry,licn < nzlicn<nz.
  ifail = 4ifail=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 11 to n.
  ifail = 5ifail=5
The input matrix is incompatible with the matrix factorized by the previous call of nag_matop_real_gen_sparse_lu (f01br) (see Section [Further Comments]).
  ifail = 6ifail=6
The input matrix is numerically singular.
W ifail = 7ifail=7
A very small pivot has been detected (see Section [Parameters], eta). The factorization has been completed but is potentially unstable.
  ifail = 8ifail=8
Duplicate elements have been found in the input matrix and the factorization has been abandoned (abort = trueabort=true on entry).

Accuracy

The factorization obtained is exact for a perturbed matrix whose (i,j)(i,j)th element differs from aijaij by less than 3ερmij3ερmij where εε is the machine precision, ρρ is the growth value returned in w(1)w1 if grow = truegrow=true, and mijmij the number of Gaussian elimination operations applied to element (i,j)(i,j).
If ρ = w(1)ρ=w1 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 44 to 77 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 (grow = truegrow=true), then the time increases by between 5%5% and 10%10%. Pivot size monitoring (eta1.0eta1.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 fill-ins. Therefore we permit additional nonzeros in positions corresponding to fill-ins.
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

function nag_matop_real_gen_sparse_lu_reuse_example
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');
irn(1:15) = [int64(1);2;2;2;3;4; ...
                    4;4;5;5;5;5; ...
                    6;6;6];
icn = zeros(150,1,'int64');
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 = [int64(1);2;2;2;3;4;4;4;5;5;5;5;6;6;6];
jvect = [int64(1);2;3;4;3;1;4;5;1;4;5;6;1;2;6];
grow = true;
[a, irn, icn, ikeep, w, idisp, ifail] = ...
    nag_matop_real_gen_sparse_lu(n, nz, a, irn, icn, abort);
a(1:15) = [10; 12; -3; -1; 15; -2; 10; -1; -1; -5; 1; -1; -1; -2; 6 ];
[aOut, w, rpmin, ifail] = ...
    nag_matop_real_gen_sparse_lu_reuse(n, a, ivect, jvect, icn, ikeep, grow, idisp)
 

aOut =

   -1.0000
   -3.0000
   -1.0000
   -2.0000
   10.0000
   15.0000
   -2.0000
    6.0000
    6.0000
   36.0000
   -1.0000
    0.0278
    1.0000
   -5.0278
    1.0000
    4.9722
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
    2.0000
   12.0000
    2.0000
   -1.0000
    6.0000
    0.2500
    2.0000
   -0.5000
   -0.5000
    1.2500
    2.0000
    2.0000


w =

   51.0000
         0
    1.0000
    5.0278
         0
         0


rpmin =

   1.0000e-04


ifail =

                    0


function f01bs_example
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');
irn(1:15) = [int64(1);2;2;2;3;4; ...
                    4;4;5;5;5;5; ...
                    6;6;6];
icn = zeros(150,1,'int64');
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 = [int64(1);2;2;2;3;4;4;4;5;5;5;5;6;6;6];
jvect = [int64(1);2;3;4;3;1;4;5;1;4;5;6;1;2;6];
grow = true;
[a, irn, icn, ikeep, w, idisp, ifail] = f01br(n, nz, a, irn, icn, abort);
a(1:15) = [10; 12; -3; -1; 15; -2; 10; -1; -1; -5; 1; -1; -1; -2; 6 ];
[aOut, w, rpmin, ifail] = f01bs(n, a, ivect, jvect, icn, ikeep, grow, idisp)
 

aOut =

   -1.0000
   -3.0000
   -1.0000
   -2.0000
   10.0000
   15.0000
   -2.0000
    6.0000
    6.0000
   36.0000
   -1.0000
    0.0278
    1.0000
   -5.0278
    1.0000
    4.9722
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
    2.0000
   12.0000
    2.0000
   -1.0000
    6.0000
    0.2500
    2.0000
   -0.5000
   -0.5000
    1.2500
    2.0000
    2.0000


w =

   51.0000
         0
    1.0000
    5.0278
         0
         0


rpmin =

   1.0000e-04


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013