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_linsys_real_sparse_fac_solve (f04ax)

Purpose

nag_linsys_real_sparse_fac_solve (f04ax) calculates the approximate solution of a set of real sparse linear equations with a single right-hand side, Ax = bAx=b or ATx = bATx=b, where AA has been factorized by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs).

Syntax

[rhs, resid] = f04ax(a, icn, ikeep, rhs, mtype, idisp, 'n', n, 'licn', licn)
[rhs, resid] = nag_linsys_real_sparse_fac_solve(a, icn, ikeep, rhs, mtype, idisp, 'n', n, 'licn', licn)

Description

To solve a system of real linear equations Ax = bAx=b or ATx = bATx=b, where AA is a general sparse matrix, AA must first be factorized by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs). nag_linsys_real_sparse_fac_solve (f04ax) then computes xx by block forward or backward substitution using simple forward and backward substitution within each diagonal block.
The method is fully described in Duff (1977).
A more recent method is available through solver function nag_sparse_direct_real_gen_solve (f11mf) and factorization function 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:     a(licn) – double array
The nonzero elements in the factorization of the matrix AA, as returned by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs).
2:     icn(licn) – int64int32nag_int array
The column indices of the nonzero elements of the factorization, as returned by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs).
3:     ikeep(5 × n5×n) – int64int32nag_int array
ikeep provides, on entry, indexing information about the factorization, as returned by nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs). Used as internal workspace prior to being restored and hence is unchanged.
4:     rhs(n) – double array
n, the dimension of the array, must satisfy the constraint n0n0.
The right-hand side vector bb.
5:     mtype – int64int32nag_int scalar
mtype specifies the task to be performed.
mtype = 1mtype=1
Solve Ax = bAx=b.
mtype1mtype1
Solve ATx = bATx=b.
6:     idisp(22) – int64int32nag_int array
The values returned in idisp by nag_matop_real_gen_sparse_lu (f01br).

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the array rhs.
nn, the order of the matrix AA.
Constraint: n0n0.
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_linsys_real_sparse_fac_solve (f04ax) is called.

Input Parameters Omitted from the MATLAB Interface

w

Output Parameters

1:     rhs(n) – double array
rhs stores the solution vector xx.
2:     resid – double scalar
The value of the maximum residual, max (|bijaijxj|)max(|bi-jaijxj|), over all the unsatisfied equations, in case nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs) has been used to factorize a singular or rectangular matrix.

Error Indicators and Warnings

If an error is detected in an input parameter nag_linsys_real_sparse_fac_solve (f04ax) will act as if a soft noisy exit has been requested (see Section [Soft Fail Option] in the (essin)).

Accuracy

The accuracy of the computed solution depends on the conditioning of the original matrix. Since nag_linsys_real_sparse_fac_solve (f04ax) is always used with either nag_matop_real_gen_sparse_lu (f01br) or nag_matop_real_gen_sparse_lu_reuse (f01bs), you are recommended to set grow = truegrow=true on entry to these functions and to examine the value of w(1)w1 on exit (see nag_matop_real_gen_sparse_lu (f01br) and nag_matop_real_gen_sparse_lu_reuse (f01bs)). For a detailed error analysis see page 17 of Duff (1977).
If storage for the original matrix is available then the error can be estimated by calculating the residual
r = bAx(or ​bATx)
r=b-Ax(or ​b-ATx)
and calling nag_linsys_real_sparse_fac_solve (f04ax) again to find a correction δδ for xx by solving
Aδ = r(or ​ATδ = r).
Aδ=r(or ​ATδ=r).

Further Comments

If the factorized form contains ττ nonzeros (idisp(2) = τidisp2=τ) then the time taken is very approximately 2τ2τ times longer than the inner loop of full matrix code. Some advantage is taken of zeros in the right-hand side when solving ATx = bATx=b (mtype1mtype1).

Example

function nag_linsys_real_sparse_fac_solve_example
n = int64(6);
nz = 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];

rhs = [15;
     12;
     18;
     3;
     -6;
     0];
mtype = int64(1);
[a, irn, icn, ikeep, w, idisp, ifail] = nag_matop_real_gen_sparse_lu(n, nz, a, irn, icn, abort);
[rhsOut, resid] = nag_linsys_real_sparse_fac_solve(a, icn, ikeep, rhs, mtype, idisp)
 

rhsOut =

     3
     3
     6
     6
     3
     1


resid =

     0


function f04ax_example
n = int64(6);
nz = 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];

rhs = [15;
     12;
     18;
     3;
     -6;
     0];
mtype = int64(1);
[a, irn, icn, ikeep, w, idisp, ifail] = f01br(n, nz, a, irn, icn, abort);
[rhsOut, resid] = f04ax(a, icn, ikeep, rhs, mtype, idisp)
 

rhsOut =

     3
     3
     6
     6
     3
     1


resid =

     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