Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_linsys_real_tridiag_fac_solve (f04le)

## Purpose

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$ by $n$ tridiagonal matrix $\left(T-\lambda I\right)$ as
 $T-λI=PLU$
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-λ Ix=y, T-λ ITx=y, Ux=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:     $\mathrm{job}$int64int32nag_int scalar
Must specify the equations to be solved.
${\mathbf{job}}=1$
The equations $\left(T-\lambda I\right)x=y$ are to be solved, but diagonal elements of $U$ are not to be perturbed.
${\mathbf{job}}=-1$
The equations $\left(T-\lambda I\right)x=y$ are to be solved and, if overflow would otherwise occur, diagonal elements of $U$ are to be perturbed. See argument tol.
${\mathbf{job}}=2$
The equations ${\left(T-\lambda I\right)}^{\mathrm{T}}x=y$ are to be solved, but diagonal elements of $U$ are not to be perturbed.
${\mathbf{job}}=-2$
The equations ${\left(T-\lambda I\right)}^{\mathrm{T}}x=y$ are to be solved and, if overflow would otherwise occur, diagonal elements of $U$ are to be perturbed. See argument tol.
${\mathbf{job}}=3$
The equations $Ux=y$ are to be solved, but diagonal elements of $U$ are not to be perturbed.
${\mathbf{job}}=-3$
The equations $Ux=y$ are to be solved and, if overflow would otherwise occur, diagonal elements of $U$ are to be perturbed. See argument tol.
2:     $\mathrm{a}\left({\mathbf{n}}\right)$ – double array
The diagonal elements of $U$ as returned by nag_linsys_real_tridiag_fac_solve (f04le).
3:     $\mathrm{b}\left({\mathbf{n}}\right)$ – double array
The elements of the first superdiagonal of $U$ as returned by nag_linsys_real_tridiag_fac_solve (f04le).
4:     $\mathrm{c}\left({\mathbf{n}}\right)$ – double array
The subdiagonal elements of $L$ as returned by nag_linsys_real_tridiag_fac_solve (f04le).
5:     $\mathrm{d}\left({\mathbf{n}}\right)$ – double array
The elements of the second superdiagonal of $U$ as returned by nag_linsys_real_tridiag_fac_solve (f04le).
6:     $\mathrm{ipiv}\left({\mathbf{n}}\right)$int64int32nag_int array
Details of the matrix $P$ as returned by nag_matop_real_gen_tridiag_lu (f01le).
7:     $\mathrm{y}\left({\mathbf{n}}\right)$ – double array
The right-hand side vector $y$.
8:     $\mathrm{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 $\epsilon ‖U‖$, where $\epsilon$ is the machine precision, but if tol is supplied as non-positive, then it is reset to $\epsilon \mathrm{max}\left|{u}_{ij}\right|$.

### Optional Input Parameters

1:     $\mathrm{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: ${\mathbf{n}}\ge 1$.

### Output Parameters

1:     $\mathrm{y}\left({\mathbf{n}}\right)$ – double array
y stores the solution vector $x$.
2:     $\mathrm{tol}$ – double scalar
If on entry tol is non-positive, it is reset as just described. Otherwise tol is unchanged.
3:     $\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:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{n}}<1$, or ${\mathbf{job}}=0$, or ${\mathbf{job}}<-3$ or ${\mathbf{job}}>3$.
${\mathbf{ifail}}>1$
Overflow would occur when computing the (${\mathbf{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.
${\mathbf{ifail}}=-99$
${\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 computed solution of the equations $\left(T-\lambda I\right)x=y$, say $\stackrel{-}{x}$, will satisfy an equation of the form
 $T-λI+Ex-=y,$
where $E$ can be expected to satisfy a bound of the form
 $E≤αεT-λI,$
$\alpha$ being a modest constant and $\epsilon$ being the machine precision. The computed solution of the equations ${\left(T-\lambda I\right)}^{\mathrm{T}}x=y$ and $Ux=y$ will satisfy similar results. The above result implies that the relative error in $\stackrel{-}{x}$ satisfies
 $x--x x- ≤cT-λIαε,$
where $c\left(T-\lambda I\right)$ is the condition number of $\left(T-\lambda I\right)$ with respect to inversion. Thus if $\left(T-\lambda I\right)$ is nearly singular, $\stackrel{-}{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.

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
 $Tx=y and TTx=y$
where
 $T= 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 and y = 2.7 -0.5 2.6 0.6 2.7 .$
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 $Tx=y$ and the second to solve ${T}^{\mathrm{T}}x=y$.
```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

```