NAG Library Routine Document

F01LEF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

+− 9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

1 Purpose

F01LEF computes an

L U

factorization of a real tridiagonal matrix, using Gaussian elimination with partial pivoting.

2 Specification

SUBROUTINE F01LEF (

N, A, LAMBDA, B, C, TOL, D, IPIV, IFAIL)

INTEGER	N, IPIV(N), IFAIL
REAL (KIND=nag_wp)	A(N), LAMBDA, B(N), C(N), TOL, D(N)

3 Description

The matrix

T - λ I

, where

T

is a real

n

n

tridiagonal matrix, is factorized as

T - λ I = P L U,

where

P

is a permutation matrix,

L

is a unit lower triangular matrix with at most one nonzero subdiagonal element per column, and

U

is an upper triangular matrix with at most two nonzero superdiagonal elements per column.

The factorization is obtained by Gaussian elimination with partial pivoting and implicit row scaling.

An indication of whether or not the matrix

T - λ I

is nearly singular is returned in the

n

th element of the array IPIV. If it is important that

T - λ I

is nonsingular, as is usually the case when solving a system of tridiagonal equations, then it is strongly recommended that

IPIV (n)

is inspected on return from F01LEF. (See the parameter IPIV and Section 8 for further details.)

The parameter

λ

is included in the routine so that F01LEF may be used, in conjunction with F04LEF, to obtain eigenvectors of

T

by inverse iteration.

4 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

5 Parameters

1: N – INTEGERInput: On entry: $n$ , the order of the matrix $T$ .
Constraint: $N \geq 1$ .
2: A(N) – REAL (KIND=nag_wp) arrayInput/Output: On entry: the diagonal elements of $T$ .

On exit: the diagonal elements of the upper triangular matrix $U$ .
3: LAMBDA – REAL (KIND=nag_wp)Input: On entry: the scalar $λ$ . F01LEF factorizes $T - λ I$ .
4: B(N) – REAL (KIND=nag_wp) arrayInput/Output: On entry: the superdiagonal elements of $T$ , stored in $B (2)$ to $B (n)$ ; $B (1)$ is not used.

On exit: the elements of the first superdiagonal of $U$ , stored in $B (2)$ to $B (n)$ .
5: C(N) – REAL (KIND=nag_wp) arrayInput/Output: On entry: the subdiagonal elements of $T$ , stored in $C (2)$ to $C (n)$ ; $C (1)$ is not used.

On exit: the subdiagonal elements of $L$ , stored in $C (2)$ to $C (n)$ .
6: TOL – REAL (KIND=nag_wp)Input: On entry: a relative tolerance used to indicate whether or not the matrix ( $T - λ I$ ) is nearly singular. TOL should normally be chosen as approximately the largest relative error in the elements of $T$ . For example, if the elements of $T$ are correct to about $4$ significant figures, then TOL should be set to about $5 \times 10^{- 4}$ . See Section 8 for further details on how TOL is used. If TOL is supplied as less than $ε$ , where $ε$ is the machine precision, then the value $ε$ is used in place of TOL.
7: D(N) – REAL (KIND=nag_wp) arrayOutput: On exit: the elements of the second superdiagonal of $U$ , stored in $D (3)$ to $D (n)$ ; $D (1)$ and $D (2)$ are not used.
8: IPIV(N) – INTEGER arrayOutput: On exit: details of the permutation matrix $P$ . If an interchange occurred at the $k$ th step of the elimination, then $IPIV (k) = 1$ , otherwise $IPIV (k) = 0$ . If a diagonal element of $U$ is small, indicating that $(T - λ I)$ is nearly singular, then the element $IPIV (n)$ is returned as positive. Otherwise $IPIV (n)$ is returned as $0$ . See Section 8 for further details. If the application is such that it is important that $(T - λ I)$ is not nearly singular, then it is strongly recommended that $IPIV (n)$ is inspected on return.
9: IFAIL – INTEGERInput/Output: On entry: IFAIL must be set to $0$ , $- 1 or 1$ . If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $- 1 or 1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is $0$ . When the value $- 1 or 1$ is used it is essential to test the value of IFAIL on exit.

On exit: $IFAIL = 0$ unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

IFAIL = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by X04AAF).

Errors or warnings detected by the routine:

$IFAIL = 1$

On entry,

N < 1

7 Accuracy

The computed factorization will satisfy the equation

P L U = (T - λ I) + E,

where

{‖E‖}_{1} \leq 9 \times \max_{i \geq j} (|l_{i j}|, {|l_{i j}|}^{2}) ε {‖T - λ I‖}_{1}

where

ε

is the machine precision.

8 Further Comments

The time taken by F01LEF is approximately proportional to

n

The factorization of a tridiagonal matrix proceeds in

(n - 1)

steps, each step eliminating one subdiagonal element of the tridiagonal matrix. In order to avoid small pivot elements and to prevent growth in the size of the elements of

L

, rows

k

and (

k + 1

) will, if necessary, be interchanged at the

k

th step prior to the elimination.

The element

IPIV (n)

returns the smallest integer,

j

, for which

|u_{j j}| \leq {‖{(T - λ I)}_{j}‖}_{1} \times TOL,

where

{‖{(T - λ I)}_{j}‖}_{1}

denotes the sum of the absolute values of the

j

th row of the matrix (

T - λ I

). If no such

j

exists, then

IPIV (n)

is returned as zero. If such a

j

exists, then

|u_{j j}|

is small and hence (

T - λ I

) is singular or nearly singular.

This routine may be followed by F04LEF to solve systems of tridiagonal equations. If you wish to solve single systems of tridiagonal equations you should be aware of F07CAF (DGTSV), which solves tridiagonal systems with a single call. F07CAF (DGTSV) requires less storage and will generally be faster than the combination of F01LEF and F04LEF, but no test for near singularity is included in F07CAF (DGTSV) and so it should only be used when the equations are known to be nonsingular.

9 Example

This example factorizes the tridiagonal matrix

T

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 then to solve the equations

T x = y

, where

y = (\begin{array}{r} 2.7 \\ - 0.5 \\ 2.6 \\ 0.6 \\ 2.7 \end{array})

by a call to F04LEF. The example program sets

TOL = 5 \times 10^{- 5}

and, of course, sets

LAMBDA = 0

NAG Library Routine DocumentF01LEF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine Document

F01LEF