Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_matop_real_symm_posdef_geneig (f01bv)

## Purpose

nag_matop_real_symm_posdef_geneig (f01bv) transforms the generalized symmetric-definite eigenproblem $Ax=\lambda {\mathbf{b}}x$ to the equivalent standard eigenproblem $Cy=\lambda y$, where $A$, ${\mathbf{b}}$ and $C$ are symmetric band matrices and ${\mathbf{b}}$ is positive definite. ${\mathbf{b}}$ must have been decomposed by nag_matop_real_symm_posdef_fac (f01bu).

## Syntax

[a, b, ifail] = f01bv(k, a, b, 'n', n, 'ma1', ma1, 'mb1', mb1)
[a, b, ifail] = nag_matop_real_symm_posdef_geneig(k, a, b, 'n', n, 'ma1', ma1, 'mb1', mb1)
Note: the interface to this routine has changed since earlier releases of the toolbox:
 At Mark 22: ma1 and mb1 were made optional

## 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 nag_matop_real_symm_posdef_fac (f01bu) as $ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}$. nag_matop_real_symm_posdef_geneig (f01bv) applies $U$, $L$ and $D$ to $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 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:     $\mathrm{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: ${\mathbf{mb1}}-1\le {\mathbf{k}}\le {\mathbf{n}}$.
2:     $\mathrm{a}\left(\mathit{lda},{\mathbf{n}}\right)$ – double array
lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge {\mathbf{ma1}}$.
The upper triangle of the $n$ by $n$ symmetric band matrix $A$, with the diagonal of the matrix stored in the $\left({m}_{A}+1\right)$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
 $* * a13 a24 a35 a46 * a12 a23 a34 a45 a56 a11 a22 a33 a44 a55 a66$
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:     $\mathrm{b}\left(\mathit{ldb},{\mathbf{n}}\right)$ – double array
ldb, the first dimension of the array, must satisfy the constraint $\mathit{ldb}\ge {\mathbf{mb1}}$.
The elements of the decomposition of matrix $B$ as returned by nag_matop_real_symm_posdef_fac (f01bu).

### Optional Input Parameters

1:     $\mathrm{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:     $\mathrm{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 ${\mathbf{ma1}}\ll {\mathbf{n}}$.
3:     $\mathrm{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: ${\mathbf{mb1}}\le {\mathbf{ma1}}$.

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},{\mathbf{n}}\right)$ – double array
Stores the corresponding elements of $C$.
2:     $\mathrm{b}\left(\mathit{ldb},{\mathbf{n}}\right)$ – double array
The elements of b will have been permuted.
3:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{mb1}}>{\mathbf{ma1}}$.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

In general the computed system is exactly congruent to a problem $\left(A+E\right)x=\lambda \left(B+F\right)x$, where $‖E‖$ and $‖F‖$ are of the order of $\epsilon \kappa \left(B\right)‖A‖$ and $\epsilon \kappa \left(B\right)‖B‖$ respectively, where $\kappa \left(B\right)$ is the condition number of $B$ with respect to inversion and $\epsilon$ 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.

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=3n/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 $Py=\lambda y$ where $P={L}^{-1}A{L}^{-\mathrm{T}}$ and $B=L{L}^{\mathrm{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}^{\mathrm{T}}$, the generalized symmetric problem reduces to the standard form
 $S-1AS-TSTx=λSTx$
and there does exist a factorization such that ${S}^{-1}A{S}^{-\mathrm{T}}$ is still of band form (see Crawford (1973)). Writing
 $C=S-1AS-T and y=STx$
the standard form is $Cy=\lambda 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 $ULD{L}^{\mathrm{T}}{U}^{\mathrm{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 $Ax=\lambda Bx$, where
 $A= 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$
 $B= 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 .$
```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

```