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)

▸▿ Contents

1 Purpose

2 Syntax

3 Description

4 References

▸▿ 5 Parameters

5.1 Compulsory Input Parameters

5.2 Optional Input Parameters

5.3 Output Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

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,

A x = b

A^{T} x = b

, where

A

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

A x = b

A^{T} x = b

, where

A

is a general sparse matrix,

A

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

x

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

A

, 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 \times 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

The right-hand side vector

b

5: $mtype$ – int64int32nag_int scalar

mtype specifies the task to be performed.

$mtype = 1$: Solve $A x = b$ .
$mtype \neq 1$: Solve $A^{T} x = b$ .

6: $idisp (2)$ – 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.
$n$ , the order of the matrix $A$ .

Constraint: $n \geq 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.

Output Parameters

1: $rhs (n)$ – double array: rhs stores the solution vector $x$ .
2: $resid$ – double scalar: The value of the maximum residual, $\max (|b_{i} - \sum_{j}^{} a_{i j} x_{j}|)$ , 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

None.

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 = true

on entry to these functions and to examine the value of

w (1)

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 = b - A x (or ​ b - A^{T} x)

and calling nag_linsys_real_sparse_fac_solve (f04ax) again to find a correction

δ

for

x

by solving

A δ = r (or ​ A^{T} δ = r) .

Further Comments

If the factorized form contains

τ

nonzeros (

idisp (2) = τ

) then the time taken is very approximately

2 τ

times longer than the inner loop of full matrix code. Some advantage is taken of zeros in the right-hand side when solving

A^{T} x = b

(

mtype \neq 1

Example

This example solves the set of linear equations

A x = b

where

A = (\begin{array}{r} 5 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & - 1 & 2 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 & 0 \\ - 2 & 0 & 0 & 1 & 1 & 0 \\ - 1 & 0 & 0 & - 1 & 2 & - 3 \\ - 1 & - 1 & 0 & 0 & 0 & 6 \end{array}) and b = (\begin{array}{r} 15 \\ 12 \\ 18 \\ 3 \\ - 6 \\ 0 \end{array}) .

The nonzero elements of

A

and indexing information are read in by the program, as described in the document for nag_matop_real_gen_sparse_lu (f01br).

Open in the MATLAB editor: f04ax_example

function f04ax_example


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

% Solve sparse system Ax = b
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');
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];

b   = [15;
       12;
       18;
        3;
       -6;
        0];

% Factorize A
abort = [true;     true;     false;     true];
[a, irn, icn, ikeep, w, idisp, ifail] = ...
f01br( ...
       n, nz, a, irn, icn, abort);

% Solve Ax = b
mtype = int64(1);
[x, resid] = f04ax( ...
                    a, icn, ikeep, b, mtype, idisp);

disp('Solution');
disp(x);

f04ax example results

Solution
     3
     3
     6
     6
     3
     1

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox