PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_matop_real_symm_posdef_geneig (f01bv)
Purpose
nag_matop_real_symm_posdef_geneig (f01bv) transforms the generalized symmetricdefinite 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 changeover point in the decomposition of
$B$ as used by
nag_matop_real_symm_posdef_fac (f01bu) and is also used as a changeover 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 bandsymmetric 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 changeover 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
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,jma1+1):j
a(ij+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$
An unexpected error has been triggered by this routine. Please
contact
NAG.
 ${\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
$\Vert E\Vert $ and
$\Vert F\Vert $ are of the order of
$\epsilon \kappa \left(B\right)\Vert A\Vert $ and
$\epsilon \kappa \left(B\right)\Vert B\Vert $ 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 wellconditioned with respect to inversion, the method, which effectively involves the inversion of
$B$, may lead to a severe loss of accuracy in wellconditioned 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=3n/4$ take $502\%$ longer.
When $B$ is positive definite and wellconditioned 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
and there does exist a factorization such that
${S}^{1}A{S}^{\mathrm{T}}$ is still of band form (see
Crawford (1973)). Writing
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 submatrix 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 righthand 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
Open in the MATLAB editor:
f01bv_example
function f01bv_example
fprintf('f01bv example results\n\n');
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];
k = int64(4);
[b, ifail] = f01bu(k, b);
[c, b, ifail] = f01bv(k, a, b);
vect = 'No vectors';
uplo = 'Upper';
q = zeros(n, n);
[~, d, e, ~, info] = f08he( ...
vect, uplo, kd, c, q);
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
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015