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_zhbgst (f08us)

## Purpose

nag_lapack_zhbgst (f08us) reduces a complex Hermitian-definite generalized eigenproblem $Az=\lambda Bz$ to the standard form $Cy=\lambda y$, where $A$ and $B$ are band matrices, $A$ is a complex Hermitian matrix, and $B$ has been factorized by nag_lapack_zpbstf (f08ut).

## Syntax

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

## Description

To reduce the complex Hermitian-definite generalized eigenproblem $Az=\lambda Bz$ to the standard form $Cy=\lambda y$, where $A$, $B$ and $C$ are banded, nag_lapack_zhbgst (f08us) must be preceded by a call to nag_lapack_zpbstf (f08ut) which computes the split Cholesky factorization of the positive definite matrix $B$: $B={S}^{\mathrm{H}}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{H}}AX$, where $X={S}^{-1}Q$ and $Q$ is a unitary 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)$ – complex 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$ Hermitian 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)$ – complex 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_zpbstf (f08ut).

### 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)$ – complex 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)$ – complex 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: rwork, 14: 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_zhbgst (f08us) 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 real floating-point operations is approximately $20{n}^{2}{k}_{B}$, when ${\mathbf{vect}}=\text{'N'}$, assuming $n\gg {k}_{A},{k}_{B}$; there are an additional $5{n}^{3}\left({k}_{B}/{k}_{A}\right)$ operations when ${\mathbf{vect}}=\text{'V'}$.
The real analogue of this function is nag_lapack_dsbgst (f08ue).

## Example

This example computes all the eigenvalues of $Az=\lambda Bz$, where
 $A = -1.13+0.00i 1.94-2.10i -1.40+0.25i 0.00+0.00i 1.94+2.10i -1.91+0.00i -0.82-0.89i -0.67+0.34i -1.40-0.25i -0.82+0.89i -1.87+0.00i -1.10-0.16i 0.00+0.00i -0.67-0.34i -1.10+0.16i 0.50+0.00i$
and
 $B = 9.89+0.00i 1.08-1.73i 0.00+0.00i 0.00+0.00i 1.08+1.73i 1.69+0.00i -0.04+0.29i 0.00+0.00i 0.00+0.00i -0.04-0.29i 2.65+0.00i -0.33+2.24i 0.00+0.00i 0.00+0.00i -0.33-2.24i 2.17+0.00i .$
Here $A$ is Hermitian, $B$ is Hermitian positive definite, and $A$ and $B$ are treated as band matrices. $B$ must first be factorized by nag_lapack_zpbstf (f08ut). The program calls nag_lapack_zhbgst (f08us) to reduce the problem to the standard form $Cy=\lambda y$, then nag_lapack_zhbtrd (f08hs) to reduce $C$ to tridiagonal form, and nag_lapack_dsterf (f08jf) to compute the eigenvalues.
```function f08us_example

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

% Sove Az = lambda Bz
% A and B are the Hermitian banded positive definite matrices:
n = 4;
ka = int64(2);
a = [ -1.13+0.00i  1.94-2.10i -1.40+0.25i  0.00+0.00i;
1.94+2.10i -1.91+0.00i -0.82-0.89i -0.67+0.34i;
-1.40-0.25i -0.82+0.89i -1.87+0.00i -1.10-0.16i;
0.00+0.00i -0.67-0.34i -1.10+0.16i  0.50+0.00i];

kb = int64(1);
b = [  9.89+0.00i  1.08-1.73i  0.00+0.00i  0.00+0.00i;
1.08+1.73i  1.69+0.00i -0.04+0.29i  0.00+0.00i;
0.00+0.00i -0.04-0.29i  2.65+0.00i -0.33+2.24i;
0.00+0.00i  0.00+0.00i -0.33-2.24i  2.17+0.00i];

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

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

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

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

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

```
```f08us example results

Eigenvalues:
-6.6089   -2.0416    0.1603    1.7712

```