NAG Library Routine Document

Integer, Intent (In)	::	n, ma1, mb1, lda, ldb, lwork
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	iwork(n)
Real (Kind=nag_wp), Intent (In)	::	relep, rmu
Real (Kind=nag_wp), Intent (Inout)	::	a(lda,n), b(ldb,n), d(30)
Real (Kind=nag_wp), Intent (Out)	::	vec(n), work(lwork)
Logical, Intent (In)	::	sym

C Header Interface

#include <nagmk26.h>

void	f02sdf_ (const Integer n, const Integer ma1, const Integer mb1, double a[], const Integer lda, double b[], const Integer ldb, const logical sym, const double relep, const double rmu, double vec[], double d[], Integer iwork[], double work[], const Integer lwork, Integer ifail)

3

Description

Given an approximation

μ

to a real eigenvalue

λ

of the generalized eigenproblem

A x = λ B x

, f02sdf attempts to compute the corresponding eigenvector by inverse iteration.

f02sdf first computes lower and upper triangular factors,

L

and

U

, of

A - μ B

, using Gaussian elimination with interchanges, and then solves the equation

U x = e

, where

e = {(1, 1, 1, \dots, 1)}^{T}

– this is the first half iteration.

There are then three possible courses of action depending on the input value of

d (1)

1.	$d (1) = 0$ . This setting should be used if $λ$ is an ill-conditioned eigenvalue (provided the matrix elements do not vary widely in order of magnitude). In this case it is essential to accept only a vector found after one half iteration, and $μ$ must be a very good approximation to $λ$ . If acceptable growth is achieved in the solution of $U x = e$ , then the normalized $x$ is accepted as the eigenvector. If not, columns of an orthogonal matrix are tried in turn in place of $e$ . If none of these give acceptable growth, the routine fails, indicating that $μ$ was not a sufficiently good approximation to $λ$ .
2.	$d (1) > 0$ . This setting should be used if $μ$ is moderately close to an eigenvalue which is not ill-conditioned (provided the matrix elements do not differ widely in order of magnitude). If acceptable growth is achieved in the solution of $U x = e$ , the normalized $x$ is accepted as the eigenvector. If not, inverse iteration is performed. Up to $30$ iterations are allowed to achieve a vector and a correction to $μ$ which together give acceptably small residuals.
3.	$d (1) < 0$ . This setting should be used if the elements of $A$ and $B$ vary widely in order of magnitude. Inverse iteration is performed, but a different convergence criterion is used.

See Section 9.3 for further details.

Note that the bandwidth of the matrix

A

must not be less than the bandwidth of

B

. If this is not so, either

A

must be filled out with zeros, or matrices

A

and

B

may be reversed and

1 / μ

supplied as an approximation to the eigenvalue

1 / λ

. Also it is assumed that

A

and

B

each have the same number of subdiagonals as superdiagonals. If this is not so, they must be filled out with zeros. If

A

and

B

are both symmetric, only the upper triangles need be supplied.

4

References

Peters G and Wilkinson J H (1979) Inverse iteration, ill-conditioned equations and Newton's method SIAM Rev. 21 339–360

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

Wilkinson J H (1972) Inverse iteration in theory and practice Symposia Mathematica Volume X 361–379 Istituto Nazionale di Alta Matematica, Monograf, Bologna

Wilkinson J H (1974) Notes on inverse iteration and ill-conditioned eigensystems Acta Univ. Carolin. Math. Phys. 1–2 173–177

Wilkinson J H (1979) Kronecker's canonical form and the

Q Z

algorithm Linear Algebra Appl. 28 285–303

5

Arguments

1: $n$ – IntegerInput

On entry:

n

, the order of the matrices

A

and

B

Constraint:

n \geq 1

2: $ma1$ – IntegerInput

On entry: the value

m_{A} + 1

, where

m_{A}

is the number of nonzero lines on each side of the diagonal of

A

. Thus the total bandwidth of

A

2 m_{A} + 1

Constraint:

1 \leq ma1 \leq n

3: $mb1$ – IntegerInput

On entry: if

mb1 \leq 0

B

is assumed to be the unit matrix. Otherwise mb1 must specify the value

m_{B} + 1

, where

m_{B}

is the number of nonzero lines on each side of the diagonal of

B

. Thus the total bandwidth of

B

2 m_{B} + 1

Constraint:

mb1 \leq ma1

4: $a (lda, n)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: the

n

n

band matrix

A

. The

m_{A}

subdiagonals must be stored in the first

m_{A}

rows of the array; the diagonal in the (

m_{A} + 1

)th row; and the

m_{A}

superdiagonals in rows

m_{A} + 2

2 m_{A} + 1

. Each row of the matrix must be stored in the corresponding column of the array. For example, if

n = 6

and

m_{A} = 2

the storage scheme is:

\begin{array}{l} * & * & a_{31} & a_{42} & a_{53} & a_{64} \\ * & a_{21} & a_{32} & a_{43} & a_{54} & a_{65} \\ a_{11} & a_{22} & a_{33} & a_{44} & a_{55} & a_{66} \\ a_{12} & a_{23} & a_{34} & a_{45} & a_{56} & * \\ a_{13} & a_{24} & a_{35} & a_{46} & * & * \end{array} .

Elements of the array marked

*

need not be set. The following code assigns the matrix elements within the band to the correct elements of the array:

   Do 20 j = 1, n
      Do 10 i = max(1,j-MA1+1), min(n,j+MA1-1)
         a(i-j+MA1,j) = matrix(j,i)
10    Continue
20 Continue

sym = .TRUE.

(i.e., both

A

and

B

are symmetric), only the lower triangle of

A

need be stored in the first ma1 rows of the array.

On exit: details of the factorization of

A - \bar{λ} B

, where

\bar{λ}

is an estimate of the eigenvalue.

5: $lda$ – IntegerInput

On entry: the first dimension of the array a as declared in the (sub)program from which f02sdf is called.

Constraint:

lda \geq 2 \times ma1 - 1

6: $b (ldb, n)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: if

mb1 > 0

, b must contain the

n

n

band matrix

B

, stored in the same way as

A

. If

sym = .TRUE.

, only the lower triangle of

B

need be stored in the first mb1 rows of the array.

mb1 \leq 0

, the array is not used.

On exit: elements in the top-left corner, and in the bottom right corner if

sym = .FALSE.

, are set to zero; otherwise the array is unchanged.

7: $ldb$ – IntegerInput

On entry: the first dimension of the array b as declared in the (sub)program from which f02sdf is called.

Constraints:

if $sym = .FALSE.$ , $ldb \geq 2 \times mb1 - 1$ ;
if $sym = .TRUE.$ , $ldb \geq mb1$ .

8: $sym$ – LogicalInput

On entry: if

sym = .TRUE.

, both

A

and

B

are assumed to be symmetric and only their upper triangles need be stored. Otherwise sym must be set to .FALSE..

9: $relep$ – Real (Kind=nag_wp)Input

On entry: the relative error of the coefficients of the given matrices

A

and

B

. If the value of relep is less than the machine precision, the machine precision is used instead.

10: $rmu$ – Real (Kind=nag_wp)Input

On entry:

μ

, an approximation to the eigenvalue for which the corresponding eigenvector is required.

11: $vec (n)$ – Real (Kind=nag_wp) arrayOutput

On exit: the eigenvector, normalized so that the largest element is unity, corresponding to the improved eigenvalue

rmu + d (30)

12: $d (30)$ – Real (Kind=nag_wp) arrayInput/Output

On entry:

d (1)

must be set to indicate the type of problem (see Section 3):

$d (1) > 0.0$: Indicates a well-conditioned eigenvalue.
$d (1) = 0.0$: Indicates an ill-conditioned eigenvalue.
$d (1) < 0.0$: Indicates that the matrices have elements varying widely in order of magnitude.

On exit: if

d (1) \neq 0.0

on entry, the successive corrections to

μ

are given in

d (i)

, for

i = 1, 2, \dots, k

, where

k + 1

is the total number of iterations performed. The final correction is also given in the last position,

d (30)

, of the array. The remaining elements of d are set to zero.

d (1) = 0.0

on entry, no corrections to

μ

are computed and

d (i)

is set to

0.0

, for

i = 1, 2, \dots, 30

. Thus in all three cases the best available approximation to the eigenvalue is

rmu + d (30)

13: $iwork (n)$ – Integer arrayWorkspace

14: $work (lwork)$ – Real (Kind=nag_wp) arrayWorkspace

15: $lwork$ – IntegerInput

On entry: the dimension of the array work as declared in the (sub)program from which f02sdf is called.

Constraints:

if $d (1) \neq 0.0$ , $lwork \geq n \times (ma1 + 1)$ ;
if $d (1) = 0.0$ , $lwork \geq 2 \times n$ .

16: $ifail$ – IntegerInput/Output

On entry: ifail must be set to

0

- 1 or 1

. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation 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 argument, 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$ ,
or	$ma1 < 1$ ,
or	$ma1 > n$ ,
or	$lda < 2 \times ma1 - 1$ ,
or	$ldb < mb1$ when $sym = .TRUE.$ ,
or	$ldb < 2 \times mb1 - 1$ when $sym = .FALSE.$ (ldb is not checked if $mb1 \leq 0$ ).

$ifail = 2$

On entry,

ma1 < mb1

. Either fill out a with zeros, or reverse the roles of a and b, and replace rmu by its reciprocal, i.e., solve

B x = λ^{- 1} A x .

$ifail = 3$

On entry,	$lwork < 2 \times n$ when $d (1) = 0.0$ ,
or	$lwork < n \times (ma1 + 1)$ when $d (1) \neq 0.0$ .

$ifail = 4$: $A$ is null. If $B$ is nonsingular, all the eigenvalues are zero and any set of n orthogonal vectors forms the eigensolution.

$ifail = 5$: $B$ is null. If $A$ is nonsingular, all the eigenvalues are infinite, and the columns of the unit matrix are eigenvectors.

$ifail = 6$

On entry,

A

and

B

are both null. The eigensolution is arbitrary.

$ifail = 7$: $d (1) \neq 0.0$ on entry and convergence is not achieved in $30$ iterations. Either the eigenvalue is ill-conditioned or rmu is a poor approximation to the eigenvalue. See Section 9.3.

$ifail = 8$: $d (1) = 0.0$ on entry and no eigenvector has been found after $\min (n, 5)$ back-substitutions. rmu is not a sufficiently good approximation to the eigenvalue.

$ifail = 9$: $d (1) < 0.0$ on entry and rmu is too inaccurate for the solution to converge.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7

Accuracy

The eigensolution is exact for some problem

(A + E) x = μ (B + F) x,

where

‖E‖, ‖F‖

are of the order of

η (‖A‖ + μ ‖B‖)

, where

η

is the value used for relep.

8

Parallelism and Performance

f02sdf is not threaded in any implementation.

9

Further Comments

9.1

Timing

The time taken by f02sdf is approximately proportional to

n {(2 m_{A} + 1)}^{2}

for factorization, and to

n (2 m_{A} + 1)

for each iteration.

9.2

Storage

The storage of the matrices

A

and

B

is designed for efficiency on a paged machine.

f02sdf will work with full matrices but it will do so inefficiently, particularly in respect of storage requirements.

9.3

Algorithmic Details

Inverse iteration is performed according to the rule

(A - μ B) y_{r + 1} = B x_{r}

x_{r + 1} = \frac{1}{α_{r + 1}} y_{r + 1}

where

α_{r + 1}

is the element of

y_{r + 1}

of largest magnitude.

Thus:

(A - μ B) x_{r + 1} = \frac{1}{α_{r + 1}} B x_{r} .

Hence the residual corresponding to

x_{r + 1}

is very small if

|α_{r + 1}|

is very large (see Peters and Wilkinson (1979)). The first half iteration,

U y_{1} = e

, corresponds to taking

L^{- 1} P B x_{0} = e

μ

is a very accurate eigenvalue, then there should always be an initial vector

x_{0}

such that one half iteration gives a small residual and thus a good eigenvector. If the eigenvalue is ill-conditioned, then second and subsequent iterated vectors may not be even remotely close to an eigenvector of a neighbouring problem (see pages 374–376 of Wilkinson (1972) and Wilkinson (1974)). In this case it is essential to accept only a vector obtained after one half iteration.

However, for well-conditioned eigenvalues, there is no loss in performing more than one iteration (see page 376 of Wilkinson (1972)), and indeed it will be necessary to iterate if

μ

is not such a good approximation to the eigenvalue. When the iteration has converged,

y_{r + 1}

will be some multiple of

x_{r}

y_{r + 1} = β_{r + 1} x_{r}

, say.

Therefore

(A - μ B) β_{r + 1} x_{r} = B x_{r},

giving

(A - (μ + \frac{1}{β_{r + 1}}) B) x_{r} = 0 .

Thus

μ + \frac{1}{β_{r + 1}}

is a better approximation to the eigenvalue.

β_{r + 1}

is obtained as the element of

y_{r + 1}

which corresponds to the element of largest magnitude,

+ 1

, in

x_{r}

. The routine terminates when

‖(A - (μ + \frac{1}{β_{r}}) B) x_{r}‖

is of the order of the machine precision relative to

‖A‖ + |μ| ‖B‖

If the elements of

A

and

B

vary widely in order of magnitude, then

‖A‖

and

‖B‖

are excessively large and a different convergence test is required. The routine terminates when the difference between successive corrections to

μ

is small relative to

μ

In practice one does not necessarily know if the given problem is well-conditioned or ill-conditioned. In order to provide some information on the condition of the eigenvalue or the accuracy of

μ

in the event of failure, successive values of

\frac{1}{β_{r}}

are stored in the vector d when

d (1)

is nonzero on input. If these values appear to be converging steadily, then it is likely that

μ

was a poor approximation to the eigenvalue and it is worth trying again with

rmu + d (30)

as the initial approximation. If the values in d vary considerably in magnitude, then the eigenvalue is ill-conditioned.

A discussion of the significance of the singularity of

A

and/or

B

is given in relation to the

Q Z

algorithm in Wilkinson (1979).

10

Example

Given the generalized eigenproblem

A x = λ B x

where

A = (\begin{array}{r} 1 & 1 & 2 \\ - 1 & 2 & 1 & 2 \\ - 1 & 3 & 1 & 2 \\ - 1 & 4 & 1 \\ - 1 & 5 \end{array}) and B = (\begin{array}{l} 5 & 1 \\ 1 & 4 & 2 \\ 2 & 3 & 2 \\ 2 & 2 & 1 \\ 1 & 1 \end{array})

find the eigenvector corresponding to the approximate eigenvalue

- 12.33

Although

B

is symmetric,

A

is not, so sym must be set to .FALSE. and all the elements of

B

in the band must be supplied to the routine.

A

(as written above) has

1

subdiagonal and

2

superdiagonals, so ma1 must be set to

3

and

A

filled out with an additional subdiagonal of zeros. Each row of the matrices is read in as data in turn.

NAG Library Routine Document

f02sdf (withdraw_real_band_geneig)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Timing

9.2 Storage

9.3 Algorithmic Details

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

9.1

Timing

9.2

Storage

9.3

Algorithmic Details

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results