as used by nag_matop_real_symm_posdef_fac (f01bu) and is also used as a change-over point in the transformations applied by this function. 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 nag_lapack_dsbtrd (f08he) and nag_lapack_dsterf (f08jf). For selected eigenvalues, use nag_lapack_dsbtrd (f08he) and nag_lapack_dstebz (f08jj).

References

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

Parameters

Compulsory Input Parameters

1: $k$ – int64int32nag_int scalar

Suggested value: the optimum value is the multiple of

m_{A}

nearest to

n / 2

k

, the change-over point in the transformations. It must be the same as the value used by nag_matop_real_symm_posdef_fac (f01bu) in the decomposition of

B

Constraint:

mb1 - 1 \leq k \leq n

2: $a (lda, n)$ – double array

lda, the first dimension of the array, must satisfy the constraint

lda \geq ma1

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:

 for j=1:n
   for i=max(1,j-ma1+1):j
     a(i-j+ma1,j) = matrix(i,j);
   end
 end

3: $b (ldb, n)$ – double array

ldb, the first dimension of the array, must satisfy the constraint

ldb \geq mb1

The elements of the decomposition of matrix

B

as returned by nag_matop_real_symm_posdef_fac (f01bu).

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the second dimension of the arrays a, b. (An error is raised if these dimensions are not equal.)
$n$ , the order of the matrices $A$ , $B$ and $C$ .
2: $ma1$ – int64int32nag_int scalar: Default: the first dimension of the array a.
$m_{A} + 1$ , where $m_{A}$ is the number of nonzero superdiagonals in $A$ . Normally $ma1 ≪ n$ .
3: $mb1$ – int64int32nag_int scalar: Default: the first dimension of the array b.
$m_{B} + 1$ , where $m_{B}$ is the number of nonzero superdiagonals in $B$ .

Constraint: $mb1 \leq ma1$ .

Output Parameters

1: $a (lda, n)$ – double array: Stores the corresponding elements of $C$ .
2: $b (ldb, n)$ – double array: The elements of b will have been permuted.
3: $ifail$ – int64int32nag_int scalar: $ifail = 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$

On entry,

mb1 > ma1

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

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.

Further Comments

The time taken by nag_matop_real_symm_posdef_geneig (f01bv) 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 function, 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}

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

Open in the MATLAB editor: f01bv_example

function f01bv_example


fprintf('f01bv example results\n\n');

% The matrices A and B
n = 9;
kd = int64(1);
a = [  0,  12,  13,  14,  15,  16,  17,  18,  19;
      11,  12,  13,  14,  15,  16,  17,  18,  19];
b = [  0,  22,  23,  24,  25,  26,  27,  28,  29;
     101, 102, 103, 104, 105, 106, 107, 108, 109];

% Decompose B
k = int64(4);
[b, ifail] = f01bu(k, b);

% Use decomposition to transform Ax = lambda Bx to Cy = lambda y
[c, b, ifail] = f01bv(k, a, b);

% Reduce C to tridiagonal form
vect = 'No vectors';
uplo = 'Upper';
q = zeros(n, n);
[~, d, e, ~, info] = f08he( ...
                            vect, uplo, kd, c, q);

% Compute eigenvalues from tridiagonal form
m1 = int64(1);
m2 = int64(3);
range = 'Indexed';
order = 'Blocked';
vl = 0;
vu = 0;
abstol = 0;
[m, nsplit, w, iblock, isplit, info] = ...
f08jj( ...
       range, order, vl, vu, m1, m2, abstol, d, e);

disp('Selected eigenvalues');
disp(w(1:m)');

f01bv example results

Selected eigenvalues
   -0.2643   -0.1530   -0.0418

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox