nag_linsys_real_tridiag_fac_solve (f04le) solves a system of tridiagonal equations following the factorization by nag_matop_real_gen_tridiag_lu (f01le). This function is intended for applications such as inverse iteration as well as straightforward linear equation applications.

Syntax

[y, tol, ifail] = f04le(job, a, b, c, d, ipiv, y, tol, 'n', n)

[y, tol, ifail] = nag_linsys_real_tridiag_fac_solve(job, a, b, c, d, ipiv, y, tol, 'n', n)

Description

Following the factorization of the

n

n

tridiagonal matrix

(T - λ I)

T - λ I = P L U

by nag_matop_real_gen_tridiag_lu (f01le), nag_linsys_real_tridiag_fac_solve (f04le) may be used to solve any of the equations

(T - λ I) x = y, {(T - λ I)}^{T} x = y, U x = y

for

x

, the choice of equation being controlled by the argument job. In each case there is an option to perturb zero or very small diagonal elements of

U

, this option being intended for use in applications such as inverse iteration.

References

Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford

Wilkinson J H and Reinsch C (1971) Handbook for Automatic Computation II, Linear Algebra Springer–Verlag

Parameters

Compulsory Input Parameters

1: $job$ – int64int32nag_int scalar

Must specify the equations to be solved.

$job = 1$: The equations $(T - λ I) x = y$ are to be solved, but diagonal elements of $U$ are not to be perturbed.
$job = - 1$: The equations $(T - λ I) x = y$ are to be solved and, if overflow would otherwise occur, diagonal elements of $U$ are to be perturbed. See argument tol.
$job = 2$: The equations ${(T - λ I)}^{T} x = y$ are to be solved, but diagonal elements of $U$ are not to be perturbed.
$job = - 2$: The equations ${(T - λ I)}^{T} x = y$ are to be solved and, if overflow would otherwise occur, diagonal elements of $U$ are to be perturbed. See argument tol.
$job = 3$: The equations $U x = y$ are to be solved, but diagonal elements of $U$ are not to be perturbed.
$job = - 3$: The equations $U x = y$ are to be solved and, if overflow would otherwise occur, diagonal elements of $U$ are to be perturbed. See argument tol.

2: $a (n)$ – double array

The diagonal elements of

U

as returned by nag_linsys_real_tridiag_fac_solve (f04le).

3: $b (n)$ – double array

The elements of the first superdiagonal of

U

as returned by nag_linsys_real_tridiag_fac_solve (f04le).

4: $c (n)$ – double array

The subdiagonal elements of

L

as returned by nag_linsys_real_tridiag_fac_solve (f04le).

5: $d (n)$ – double array

The elements of the second superdiagonal of

U

as returned by nag_linsys_real_tridiag_fac_solve (f04le).

6: $ipiv (n)$ – int64int32nag_int array

Details of the matrix

P

as returned by nag_matop_real_gen_tridiag_lu (f01le).

7: $y (n)$ – double array

The right-hand side vector

y

8: $tol$ – double scalar

The minimum perturbation to be made to very small diagonal elements of

U

. tol is only referenced when job is negative. tol should normally be chosen as about

ε ‖U‖

, where

ε

is the machine precision, but if tol is supplied as non-positive, then it is reset to

ε \max |u_{i j}|

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the dimension of the arrays a, b, c, d, ipiv, y. (An error is raised if these dimensions are not equal.)
$n$ , the order of the matrix $T$ .

Constraint: $n \geq 1$ .

Output Parameters

1: $y (n)$ – double array: y stores the solution vector $x$ .
2: $tol$ – double scalar: If on entry tol is non-positive, it is reset as just described. Otherwise tol is unchanged.
3: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

$ifail = 1$

On entry,	$n < 1$ ,
or	$job = 0$ ,
or	$job < - 3$ or $job > 3$ .

$ifail > 1$: Overflow would occur when computing the ( $ifail - 1$ )th element of the solution vector $x$ . This can only occur when job is supplied as positive and either means that a diagonal element of $U$ is very small or that elements of the right-hand side vector $y$ are very large.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

The computed solution of the equations

(T - λ I) x = y

, say

\bar{x}

, will satisfy an equation of the form

(T - λ I + E) \bar{x} = y,

where

E

can be expected to satisfy a bound of the form

‖E‖ \leq α ε ‖T - λ I‖,

α

being a modest constant and

ε

being the machine precision. The computed solution of the equations

{(T - λ I)}^{T} x = y

and

U x = y

will satisfy similar results. The above result implies that the relative error in

\bar{x}

satisfies

\frac{‖\bar{x} - x‖}{‖\bar{x}‖} \leq c (T - λ I) α ε,

where

c (T - λ I)

is the condition number of

(T - λ I)

with respect to inversion. Thus if

(T - λ I)

is nearly singular,

\bar{x}

can be expected to have a large relative error. Note that nag_matop_real_gen_tridiag_lu (f01le) incorporates a test for near singularity.

Further Comments

The time taken by nag_linsys_real_tridiag_fac_solve (f04le) is approximately proportional to

n

If you have single systems of tridiagonal equations to solve you are advised that nag_lapack_dgtsv (f07ca) requires less storage and will normally be faster than the combination of nag_matop_real_gen_tridiag_lu (f01le) and nag_linsys_real_tridiag_fac_solve (f04le), but nag_lapack_dgtsv (f07ca) does not incorporate a test for near singularity.

Example

This example solves the two sets of tridiagonal equations

T x = y and T^{T} x = y

where

T = (\begin{array}{r} 3.0 & 2.1 & 0 & 0 & 0 \\ 3.4 & 2.3 & - 1.0 & 0 & 0 \\ 0 & 3.6 & - 5.0 & 1.9 & 0 \\ 0 & 0 & 7.0 & - 0.9 & 8.0 \\ 0 & 0 & 0 & - 6.0 & 7.1 \end{array}) and y = (\begin{array}{r} 2.7 \\ - 0.5 \\ 2.6 \\ 0.6 \\ 2.7 \end{array}) .

The example program first calls nag_matop_real_gen_tridiag_lu (f01le) to factorize

T

and then two calls are made to nag_linsys_real_tridiag_fac_solve (f04le), one to solve

T x = y

and the second to solve

T^{T} x = y

Open in the MATLAB editor: f04le_example

function f04le_example


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

% Tridiagonal matrix T stored by diagonals
b =       [0     2.1    -1      1.9     8];
a =       [3     2.3    -5     -0.9     7.1];
c = [0     3.4   3.6     7     -6];

% LU factorize T
lambda = 0;
tol = 5e-05;
[D, U1, L, U2, ipiv, ifail] = f01le( ...
                                     a, lambda, b, c, tol);

% Solve Tx = y, then T'x = y where 
y = [2.7  -0.5   2.6   0.6   2.7 ];
% using factorized form of T
job = int64(1);
[x, ~, ifail] = f04le( ...
                       job, D, U1, L, U2, ipiv, y, tol);

disp('Solution vector for Tx = y:');
disp(x);

job = int64(2);
[x, ~, ifail] = f04le( ...
                       job, D, U1, L, U2, ipiv, y, tol);

disp('Solution vector for T''x = y:');
disp(x);

f04le example results

Solution vector for Tx = y:
   -4.0000    7.0000    3.0000   -4.0000   -3.0000

Solution vector for T'x = y:
   -4.6304    4.8798   -0.5554    0.6718   -0.3767

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

Chapter Contents

Chapter Introduction

NAG Toolbox