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$n$ by n$n$ tridiagonal matrix (TλI)$\left(T-\lambda I\right)$ as
 T − λI = PLU $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 − λI)x = y,   (T − λI)Tx = y,   Ux = y $(T-λ I)x=y, (T-λ I)Tx=y, Ux=y$
for x$x$, the choice of equation being controlled by the parameter job. In each case there is an option to perturb zero or very small diagonal elements of U$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${\mathbf{job}}=1$
The equations (TλI)x = y$\left(T-\lambda I\right)x=y$ are to be solved, but diagonal elements of U$U$ are not to be perturbed.
job = 1${\mathbf{job}}=-1$
The equations (TλI)x = y$\left(T-\lambda I\right)x=y$ are to be solved and, if overflow would otherwise occur, diagonal elements of U$U$ are to be perturbed. See parameter tol.
job = 2${\mathbf{job}}=2$
The equations (TλI)Tx = y${\left(T-\lambda I\right)}^{\mathrm{T}}x=y$ are to be solved, but diagonal elements of U$U$ are not to be perturbed.
job = 2${\mathbf{job}}=-2$
The equations (TλI)Tx = y${\left(T-\lambda I\right)}^{\mathrm{T}}x=y$ are to be solved and, if overflow would otherwise occur, diagonal elements of U$U$ are to be perturbed. See parameter tol.
job = 3${\mathbf{job}}=3$
The equations Ux = y$Ux=y$ are to be solved, but diagonal elements of U$U$ are not to be perturbed.
job = 3${\mathbf{job}}=-3$
The equations Ux = y$Ux=y$ are to be solved and, if overflow would otherwise occur, diagonal elements of U$U$ are to be perturbed. See parameter tol.
2:     a(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
The diagonal elements of U$U$ as returned by nag_linsys_real_tridiag_fac_solve (f04le).
3:     b(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
The elements of the first superdiagonal of U$U$ as returned by nag_linsys_real_tridiag_fac_solve (f04le).
4:     c(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
The subdiagonal elements of L$L$ as returned by nag_linsys_real_tridiag_fac_solve (f04le).
5:     d(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
The elements of the second superdiagonal of U$U$ as returned by nag_linsys_real_tridiag_fac_solve (f04le).
6:     ipiv(n) – int64int32nag_int array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
Details of the matrix P$P$ as returned by nag_matop_real_gen_tridiag_lu (f01le).
7:     y(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
The right-hand side vector y$y$.
8:     tol – double scalar
The minimum perturbation to be made to very small diagonal elements of U$U$. tol is only referenced when job is negative. tol should normally be chosen as about εU$\epsilon ‖U‖$, where ε$\epsilon$ is the machine precision, but if tol is supplied as non-positive, then it is reset to εmax|uij|$\epsilon \mathrm{max}|{u}_{ij}|$.

### 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$n$, the order of the matrix T$T$.
Constraint: n1${\mathbf{n}}\ge 1$.

None.

### Output Parameters

1:     y(n) – double array
y stores the solution vector x$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
${\mathrm{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:
ifail = 1${\mathbf{ifail}}=1$
 On entry, n < 1${\mathbf{n}}<1$, or job = 0${\mathbf{job}}=0$, or job < − 3${\mathbf{job}}<-3$ or job > 3${\mathbf{job}}>3$.
ifail > 1${\mathbf{ifail}}>1$
Overflow would occur when computing the (ifail1${\mathbf{ifail}}-1$)th element of the solution vector x$x$. This can only occur when job is supplied as positive and either means that a diagonal element of U$U$ is very small or that elements of the right-hand side vector y$y$ are very large.

## Accuracy

The computed solution of the equations (TλI)x = y$\left(T-\lambda I\right)x=y$, say x$\stackrel{-}{x}$, will satisfy an equation of the form
 (T − λI + E)x = y, $(T-λI+E)x-=y,$
where E$E$ can be expected to satisfy a bound of the form
 ‖E‖ ≤ αε‖T − λI‖, $‖E‖≤αε‖T-λI‖,$
α$\alpha$ being a modest constant and ε$\epsilon$ being the machine precision. The computed solution of the equations (TλI)Tx = y${\left(T-\lambda I\right)}^{\mathrm{T}}x=y$ and Ux = y$Ux=y$ will satisfy similar results. The above result implies that the relative error in x$\stackrel{-}{x}$ satisfies
 (‖x − x‖)/(‖x‖) ≤ c(T − λI)αε, $‖x--x‖ ‖x-‖ ≤c(T-λI)αε,$
where c(TλI)$c\left(T-\lambda I\right)$ is the condition number of (TλI)$\left(T-\lambda I\right)$ with respect to inversion. Thus if (TλI)$\left(T-\lambda I\right)$ is nearly singular, x$\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$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

```function nag_linsys_real_tridiag_fac_solve_example
job = int64(1);
a = [3;
2.3;
-5;
-0.9;
7.1];
lambda = 0;
b = [0;
2.1;
-1;
1.9;
8];
c = [0;
3.4;
3.6;
7;
-6];
tol = 5e-05;
ipiv = [int64(0);1;1;1;0];
y = [2.7;
-0.5;
2.6;
0.6;
2.7];
tol = 5e-05;
[a, b, c, d, ipiv, ifail] = nag_matop_real_gen_tridiag_lu(a, lambda, b, c, tol);
[yOut, tolOut, ifail] = nag_linsys_real_tridiag_fac_solve(job, a, b, c, d, ipiv, y, tol)
```
```

yOut =

-4.0000
7.0000
3.0000
-4.0000
-3.0000

tolOut =

5.0000e-05

ifail =

0

```
```function f04le_example
job = int64(1);
a = [3;
2.3;
-5;
-0.9;
7.1];
lambda = 0;
b = [0;
2.1;
-1;
1.9;
8];
c = [0;
3.4;
3.6;
7;
-6];
tol = 5e-05;
ipiv = [int64(0);1;1;1;0];
y = [2.7;
-0.5;
2.6;
0.6;
2.7];
tol = 5e-05;
[a, b, c, d, ipiv, ifail] = f01le(a, lambda, b, c, tol);
[yOut, tolOut, ifail] = f04le(job, a, b, c, d, ipiv, y, tol)
```
```

yOut =

-4.0000
7.0000
3.0000
-4.0000
-3.0000

tolOut =

5.0000e-05

ifail =

0

```