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_lapack_dsbgst (f08ue)

## Purpose

nag_lapack_dsbgst (f08ue) reduces a real symmetric-definite generalized eigenproblem $Az=\lambda Bz$ to the standard form $Cy=\lambda y$, where $A$ and $B$ are band matrices, $A$ is a real symmetric matrix, and $B$ has been factorized by nag_lapack_dpbstf (f08uf).

## Syntax

[ab, x, info] = f08ue(vect, uplo, ka, kb, ab, bb, 'n', n)
[ab, x, info] = nag_lapack_dsbgst(vect, uplo, ka, kb, ab, bb, 'n', n)

## Description

To reduce the real symmetric-definite generalized eigenproblem $Az=\lambda Bz$ to the standard form $Cy=\lambda y$, where $A$, $B$ and $C$ are banded, nag_lapack_dsbgst (f08ue) must be preceded by a call to nag_lapack_dpbstf (f08uf) which computes the split Cholesky factorization of the positive definite matrix $B$: $B={S}^{\mathrm{T}}S$. The split Cholesky factorization, compared with the ordinary Cholesky factorization, allows the work to be approximately halved.
This function overwrites $A$ with $C={X}^{\mathrm{T}}AX$, where $X={S}^{-1}Q$ and $Q$ is a orthogonal matrix chosen (implicitly) to preserve the bandwidth of $A$. The function also has an option to allow the accumulation of $X$, and then, if $z$ is an eigenvector of $C$, $Xz$ is an eigenvector of the original system.

## References

Crawford C R (1973) Reduction of a band-symmetric generalized eigenvalue problem Comm. ACM 16 41–44
Kaufman L (1984) Banded eigenvalue solvers on vector machines ACM Trans. Math. Software 10 73–86

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{vect}$ – string (length ≥ 1)
Indicates whether $X$ is to be returned.
${\mathbf{vect}}=\text{'N'}$
$X$ is not returned.
${\mathbf{vect}}=\text{'V'}$
$X$ is returned.
Constraint: ${\mathbf{vect}}=\text{'N'}$ or $\text{'V'}$.
2:     $\mathrm{uplo}$ – string (length ≥ 1)
Indicates whether the upper or lower triangular part of $A$ is stored.
${\mathbf{uplo}}=\text{'U'}$
The upper triangular part of $A$ is stored.
${\mathbf{uplo}}=\text{'L'}$
The lower triangular part of $A$ is stored.
Constraint: ${\mathbf{uplo}}=\text{'U'}$ or $\text{'L'}$.
3:     $\mathrm{ka}$int64int32nag_int scalar
If ${\mathbf{uplo}}=\text{'U'}$, the number of superdiagonals, ${k}_{a}$, of the matrix $A$.
If ${\mathbf{uplo}}=\text{'L'}$, the number of subdiagonals, ${k}_{a}$, of the matrix $A$.
Constraint: ${\mathbf{ka}}\ge 0$.
4:     $\mathrm{kb}$int64int32nag_int scalar
If ${\mathbf{uplo}}=\text{'U'}$, the number of superdiagonals, ${k}_{b}$, of the matrix $B$.
If ${\mathbf{uplo}}=\text{'L'}$, the number of subdiagonals, ${k}_{b}$, of the matrix $B$.
Constraint: ${\mathbf{ka}}\ge {\mathbf{kb}}\ge 0$.
5:     $\mathrm{ab}\left(\mathit{ldab},:\right)$ – double array
The first dimension of the array ab must be at least ${\mathbf{ka}}+1$.
The second dimension of the array ab must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The upper or lower triangle of the $n$ by $n$ symmetric band matrix $A$.
The matrix is stored in rows $1$ to ${k}_{a}+1$, more precisely,
• if ${\mathbf{uplo}}=\text{'U'}$, the elements of the upper triangle of $A$ within the band must be stored with element ${A}_{ij}$ in ${\mathbf{ab}}\left({k}_{a}+1+i-j,j\right)\text{​ for ​}\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,j-{k}_{a}\right)\le i\le j$;
• if ${\mathbf{uplo}}=\text{'L'}$, the elements of the lower triangle of $A$ within the band must be stored with element ${A}_{ij}$ in ${\mathbf{ab}}\left(1+i-j,j\right)\text{​ for ​}j\le i\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(n,j+{k}_{a}\right)\text{.}$
6:     $\mathrm{bb}\left(\mathit{ldbb},:\right)$ – double array
The first dimension of the array bb must be at least ${\mathbf{kb}}+1$.
The second dimension of the array bb must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The banded split Cholesky factor of $B$ as specified by uplo, n and kb and returned by nag_lapack_dpbstf (f08uf).

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the second dimension of the arrays ab, bb.
$n$, the order of the matrices $A$ and $B$.
Constraint: ${\mathbf{n}}\ge 0$.

### Output Parameters

1:     $\mathrm{ab}\left(\mathit{ldab},:\right)$ – double array
The first dimension of the array ab will be ${\mathbf{ka}}+1$.
The second dimension of the array ab will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The upper or lower triangle of ab stores the corresponding upper or lower triangle of $C$ as specified by uplo.
2:     $\mathrm{x}\left(\mathit{ldx},:\right)$ – double array
The first dimension, $\mathit{ldx}$, of the array x will be
• if ${\mathbf{vect}}=\text{'V'}$, $\mathit{ldx}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$;
• if ${\mathbf{vect}}=\text{'N'}$, $\mathit{ldx}=1$.
The second dimension of the array x will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$ if ${\mathbf{vect}}=\text{'V'}$ and at least $1$ if ${\mathbf{vect}}=\text{'N'}$.
The $n$ by $n$ matrix $X={S}^{-1}Q$, if ${\mathbf{vect}}=\text{'V'}$.
If ${\mathbf{vect}}=\text{'N'}$, x is not referenced.
3:     $\mathrm{info}$int64int32nag_int scalar
${\mathbf{info}}=0$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

${\mathbf{info}}=-i$
If ${\mathbf{info}}=-i$, parameter $i$ had an illegal value on entry. The parameters are numbered as follows:
1: vect, 2: uplo, 3: n, 4: ka, 5: kb, 6: ab, 7: ldab, 8: bb, 9: ldbb, 10: x, 11: ldx, 12: work, 13: info.
It is possible that info refers to a parameter that is omitted from the MATLAB interface. This usually indicates that an error in one of the other input parameters has caused an incorrect value to be inferred.

## Accuracy

Forming the reduced matrix $C$ is a stable procedure. However it involves implicit multiplication by ${B}^{-1}$. When nag_lapack_dsbgst (f08ue) is used as a step in the computation of eigenvalues and eigenvectors of the original problem, there may be a significant loss of accuracy if $B$ is ill-conditioned with respect to inversion.

The total number of floating-point operations is approximately $6{n}^{2}{k}_{B}$, when ${\mathbf{vect}}=\text{'N'}$, assuming $n\gg {k}_{A},{k}_{B}$; there are an additional $\left(3/2\right){n}^{3}\left({k}_{B}/{k}_{A}\right)$ operations when ${\mathbf{vect}}=\text{'V'}$.
The complex analogue of this function is nag_lapack_zhbgst (f08us).

## Example

This example computes all the eigenvalues of $Az=\lambda Bz$, where
 $A = 0.24 0.39 0.42 0.00 0.39 -0.11 0.79 0.63 0.42 0.79 -0.25 0.48 0.00 0.63 0.48 -0.03 and B= 2.07 0.95 0.00 0.00 0.95 1.69 -0.29 0.00 0.00 -0.29 0.65 -0.33 0.00 0.00 -0.33 1.17 .$
Here $A$ is symmetric, $B$ is symmetric positive definite, and $A$ and $B$ are treated as band matrices. $B$ must first be factorized by nag_lapack_dpbstf (f08uf). The program calls nag_lapack_dsbgst (f08ue) to reduce the problem to the standard form $Cy=\lambda y$, then nag_lapack_dsbtrd (f08he) to reduce $C$ to tridiagonal form, and nag_lapack_dsterf (f08jf) to compute the eigenvalues.
```function f08ue_example

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

% Sove Az = lambda Bz
% A and B are the symmetric banded positive definite matrices:
n = 4;
% A has 2 off-diagonals
ka = int64(2);
a = [ 0.24,  0.39,  0.42,  0.00;
0.39, -0.11,  0.79,  0.63;
0.42,  0.79, -0.25,  0.48;
0.00,  0.63,  0.48, -0.03];
% B has 1 off-diagonal
kb = int64(1);
b = [ 2.07   0.95   0.00   0.00;
0.95   1.69  -0.29   0.00;
0.00  -0.29   0.65  -0.33;
0.00   0.00  -0.33   1.17];

% Convert to general banded format ...
[~, ab, ifail] = f01zc( ...
'P', ka, ka, a, zeros(ka+ka+1,n));
[~, bb, ifail] = f01zc( ...
'P', kb, kb, b, zeros(kb+kb+1,n));
% ... and chop to give 'Upper' symmetric banded format
ab = ab(1:ka+1,1:n);
bb = bb(1:kb+1,1:n);

% Factorize B
uplo = 'Upper';
[ub, info] = f08uf( ...
uplo, kb, bb);

% Reduce problem to standard form Cy = lambda*y
vect = 'N';
[cb, x, info] = f08ue( ...
vect, uplo, ka, kb, ab, ub);

% Find eigenvalues lambda
jobz = 'No Vectors';
[~, w, ~, info] = f08ha( ...
jobz, uplo, ka, cb);

disp('Eigenvalues:');
disp(w');

```
```f08ue example results

Eigenvalues:
-0.8305   -0.6401    0.0992    1.8525

```