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 \times 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 matrix elements within the band can be assigned to the correct elements of the array using the following code:

   Do j = 1, n
      Do i = max(1,j-ma1+1), J
         a(i-j+ma1,j) = matrix (i,j)
      End Do
   End Do

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

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

−1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

−1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

0

is recommended. 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

Background information to multithreading can be found in the Multithreading documentation.

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}) .

f01bv: FL CL CPP AD PY MB

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)