Integer, Intent (In)	::	n, ma1, mb1, m3, k, lda, ldb, ldv
Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), Intent (Inout)	::	a(lda,n), b(ldb,n), v(ldv,m3)
Real (Kind=nag_wp), Intent (Out)	::	w(m3)

C Header Interface

#include <nag.h>

void	f01bvf_ (const Integer n, const Integer ma1, const Integer mb1, const Integer m3, const Integer k, double a[], const Integer lda, double b[], const Integer ldb, double v[], const Integer ldv, double w[], Integer *ifail)

The routine may be called by the names f01bvf or nagf_matop_real_symm_posdef_geneig.

3 Description

A

is a symmetric band matrix of order

n

and bandwidth

2 m_{A} + 1

. The positive definite symmetric band matrix

B

, of order

n

and bandwidth

2 m_{B} + 1

, must have been previously decomposed by f01buf as

U L D L^{T} U^{T}

. f01bvf applies

U

L

and

D

A

m_{A}

rows at a time, restoring the band form of

A

at each stage by plane rotations. The argument

k

defines the change-over point in the decomposition of

B

as used by f01buf and is also used as a change-over point in the transformations applied by this routine. For maximum efficiency,

k

should be chosen to be the multiple of

m_{A}

nearest to

n / 2

. The resulting symmetric band matrix

C

is overwritten on a. The eigenvalues of

C

, and thus of the original problem, may be found using f08hef and f08jff. For selected eigenvalues, use f08hef and f08jjf.

4 References

Crawford C R (1973) Reduction of a band-symmetric generalized eigenvalue problem Comm. ACM 16 41–44

5 Arguments

1: $n$ – Integer Input

On entry:

n

, the order of the matrices

A

B

and

C

2: $ma1$ – Integer Input

On entry:

m_{A} + 1

, where

m_{A}

is the number of nonzero superdiagonals in

A

. Normally

ma1 ≪ n

3: $mb1$ – Integer Input

On entry:

m_{B} + 1

, where

m_{B}

is the number of nonzero superdiagonals in

B

Constraint:

mb1 \leq ma1

4: $m3$ – Integer Input

On entry: the value of

3 m_{A} + m_{B}

5: $k$ – Integer Input

On entry:

k

, the change-over point in the transformations. It must be the same as the value used by f01buf in the decomposition of

B

Suggested value: the optimum value is the multiple of

m_{A}

nearest to

n / 2

Constraint:

mb1 - 1 \leq k \leq n

6: $a (lda, n)$ – Real (Kind=nag_wp) array Input/Output

On entry: the upper triangle of the

n

n

symmetric band matrix

A

, with the diagonal of the matrix stored in the

(m_{A} + 1)

th row of the array, and the

m_{A}

superdiagonals within the band stored in the first

m_{A}

rows of the array. Each column of the matrix is stored in the corresponding column of the array. For example, if

n = 6

and

m_{A} = 2

, the storage scheme is

\begin{matrix} * & * & a_{13} & a_{24} & a_{35} & a_{46} \\ * & a_{12} & a_{23} & a_{34} & a_{45} & a_{56} \\ a_{11} & a_{22} & a_{33} & a_{44} & a_{55} & a_{66} \end{matrix}

Elements in the top left corner of the array 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), J
         A(I-J+MA1,J) = matrix (I,J)
10    CONTINUE
20 CONTINUE

On exit: is overwritten by the corresponding elements of

C

7: $lda$ – Integer Input

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

Constraint:

lda \geq ma1

8: $b (ldb, n)$ – Real (Kind=nag_wp) array Input/Output

On entry: the elements of the decomposition of matrix

B

as returned by f01buf.

On exit: the elements of b will have been permuted.

9: $ldb$ – Integer Input

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

Constraint:

ldb \geq mb1

10: $v (ldv, m3)$ – Real (Kind=nag_wp) array Workspace

11: $ldv$ – Integer Input

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

Constraint:

ldv \geq m_{A} + m_{B}

12: $w (m3)$ – Real (Kind=nag_wp) array Workspace

13: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

- 1 or 1

. If you are unfamiliar with this argument you should refer to Section 4 in the Introduction to the NAG Library FL Interface 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, $mb1 = 〈value〉$ and $ma1 = 〈value〉$ .
Constraint: $mb1 \leq ma1$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

In general the computed system is exactly congruent to a problem

(A + E) x = λ (B + F) x

, where

‖E‖

and

‖F‖

are of the order of

ε κ (B) ‖A‖

and

ε κ (B) ‖B‖

respectively, where

κ (B)

is the condition number of

B

with respect to inversion and

ε

is the machine precision. This means that when

B

is positive definite but not well-conditioned with respect to inversion, the method, which effectively involves the inversion of

B

, may lead to a severe loss of accuracy in well-conditioned eigenvalues.

8 Parallelism and Performance

f01bvf is not threaded in any implementation.

9 Further Comments

The time taken by f01bvf is approximately proportional to

n^{2} m_{B}^{2}

and the distance of

k

from

n / 2

, e.g.,

k = n / 4

and

k = 3 n / 4

take

502 %

longer.

When

B

is positive definite and well-conditioned with respect to inversion, the generalized symmetric eigenproblem can be reduced to the standard symmetric problem

P y = λ y

where

P = L^{- 1} A L^{- T}

and

B = L L^{T}

, the Cholesky factorization.

When

A

and

B

are of band form, especially if the bandwidth is small compared with the order of the matrices, storage considerations may rule out the possibility of working with

P

since it will be a full matrix in general. However, for any factorization of the form

B = S S^{T}

, the generalized symmetric problem reduces to the standard form

S^{- 1} A S^{- T} (S^{T} x) = λ (S^{T} x)

and there does exist a factorization such that

S^{- 1} A S^{- T}

is still of band form (see Crawford (1973)). Writing

C = S^{- 1} A S^{- T} and y = S^{T} x

the standard form is

C y = λ y

and the bandwidth of

C

is the maximum bandwidth of

A

and

B

Each stage in the transformation consists of two phases. The first reduces a leading principal sub-matrix of

B

to the identity matrix and this introduces nonzero elements outside the band of

A

. In the second, further transformations are applied which leave the reduced part of

B

unaltered and drive the extra elements upwards and off the top left corner of

A

. Alternatively,

B

may be reduced to the identity matrix starting at the bottom right-hand corner and the extra elements introduced in

A

can be driven downwards.

The advantage of the

U L D L^{T} U^{T}

decomposition of

B

is that no extra elements have to be pushed over the whole length of

A

. If

k

is taken as approximately

n / 2

, the shifting is limited to halfway. At each stage the size of the triangular bumps produced in

A

depends on the number of rows and columns of

B

which are eliminated in the first phase and on the bandwidth of

B

. The number of rows and columns over which these triangles are moved at each step in the second phase is equal to the bandwidth of

A

In this routine, a is defined as being at least as wide as

B

and must be filled out with zeros if necessary as it is overwritten with

C

. The number of rows and columns of

B

which are effectively eliminated at each stage is

m_{A}

10 Example

This example finds the three smallest eigenvalues of

A x = λ B x

, where

A = (\begin{array}{r} 11 & 12 \\ 12 & 12 & 13 \\ 13 & 13 & 14 \\ 14 & 14 & 15 \\ 15 & 15 & 16 \\ 16 & 16 & 17 \\ 17 & 17 & 18 \\ 18 & 18 & 19 \\ 19 & 19 \end{array})

B = (\begin{array}{r} 101 & 22 \\ 22 & 102 & 23 \\ 23 & 103 & 24 \\ 24 & 104 & 25 \\ 25 & 105 & 26 \\ 26 & 106 & 27 \\ 27 & 107 & 28 \\ 28 & 108 & 29 \\ 29 & 109 \end{array}) .

Interfaces: FL CL AD

NAG FL Interface Introduction

F01 (Matop) Chapter Contents

F01 (Matop) Chapter Introduction

f01bv: FL CL AD

NAG FL Interfacef01bvf (real_​symm_​posdef_​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

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
f01bvf (real_symm_posdef_geneig)