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 = λbx$Ax=\lambda {\mathbf{b}}x$ to the equivalent standard eigenproblem Cy = λy$Cy=\lambda y$, where A$A$, b${\mathbf{b}}$ and C$C$ are symmetric band matrices and b${\mathbf{b}}$ is positive definite. b${\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:
Mark 22: ma1, mb1 have been made optional
.

## Description

A$A$ is a symmetric band matrix of order n$n$ and bandwidth 2mA + 1$2{m}_{A}+1$. The positive definite symmetric band matrix B$B$, of order n$n$ and bandwidth 2mB + 1$2{m}_{B}+1$, must have been previously decomposed by nag_matop_real_symm_posdef_fac (f01bu) as ULDLTUT$ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}$. nag_matop_real_symm_posdef_geneig (f01bv) applies U$U$, L$L$ and D$D$ to A$A$, mA${m}_{A}$ rows at a time, restoring the band form of A$A$ at each stage by plane rotations. The parameter k$k$ defines the change-over point in the decomposition of B$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$k$ should be chosen to be the multiple of mA${m}_{A}$ nearest to n / 2$n/2$. The resulting symmetric band matrix C$C$ is overwritten on a. The eigenvalues of C$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
k$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$B$.
Constraint: mb11kn${\mathbf{mb1}}-1\le {\mathbf{k}}\le {\mathbf{n}}$.
2:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldama1$\mathit{lda}\ge {\mathbf{ma1}}$.
The upper triangle of the n$n$ by n$n$ symmetric band matrix A$A$, with the diagonal of the matrix stored in the (mA + 1)$\left({m}_{A}+1\right)$th row of the array, and the mA${m}_{A}$ superdiagonals within the band stored in the first mA${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$n=6$ and mA = 2${m}_{A}=2$, the storage scheme is
 * * a13 a24 a35 a46 * a12 a23 a34 a45 a56 a11 a22 a33 a44 a55 a66
$* * 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:     b(ldb,n) – double array
ldb, the first dimension of the array, must satisfy the constraint ldbmb1$\mathit{ldb}\ge {\mathbf{mb1}}$.
The elements of the decomposition of matrix B$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$n$, the order of the matrices A$A$, B$B$ and C$C$.
2:     ma1 – int64int32nag_int scalar
Default: The first dimension of the array a.
mA + 1${m}_{A}+1$, where mA${m}_{A}$ is the number of nonzero superdiagonals in A$A$. Normally ma1n${\mathbf{ma1}}\ll {\mathbf{n}}$.
3:     mb1 – int64int32nag_int scalar
Default: The first dimension of the array b.
mB + 1${m}_{B}+1$, where mB${m}_{B}$ is the number of nonzero superdiagonals in B$B$.
Constraint: ${\mathbf{mb1}}\le {\mathbf{ma1}}$.

### Input Parameters Omitted from the MATLAB Interface

m3 lda ldb v ldv w

### Output Parameters

1:     a(lda,n) – double array
ldama1$\mathit{lda}\ge {\mathbf{ma1}}$.
Stores the corresponding elements of C$C$.
2:     b(ldb,n) – double array
ldbmb1$\mathit{ldb}\ge {\mathbf{mb1}}$.
The elements of b will have been permuted.
3:     ifail – int64int32nag_int scalar
${\mathrm{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:
ifail = 1${\mathbf{ifail}}=1$
 On entry, ${\mathbf{mb1}}>{\mathbf{ma1}}$.

## Accuracy

In general the computed system is exactly congruent to a problem (A + E)x = λ(B + F)x$\left(A+E\right)x=\lambda \left(B+F\right)x$, where E$‖E‖$ and F$‖F‖$ are of the order of εκ(B)A$\epsilon \kappa \left(B\right)‖A‖$ and εκ(B)B$\epsilon \kappa \left(B\right)‖B‖$ respectively, where κ(B)$\kappa \left(B\right)$ is the condition number of B$B$ with respect to inversion and ε$\epsilon$ is the machine precision. This means that when B$B$ is positive definite but not well-conditioned with respect to inversion, the method, which effectively involves the inversion of B$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 n2mB2${n}^{2}{m}_{B}^{2}$ and the distance of k$k$ from n / 2$n/2$, e.g., k = n / 4$k=n/4$ and k = 3n / 4$k=3n/4$ take 502%$502%$ longer.
When B$B$ is positive definite and well-conditioned with respect to inversion, the generalized symmetric eigenproblem can be reduced to the standard symmetric problem Py = λy$Py=\lambda y$ where P = L1ALT$P={L}^{-1}A{L}^{-\mathrm{T}}$ and B = LLT$B=L{L}^{\mathrm{T}}$, the Cholesky factorization.
When A$A$ and B$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$P$ since it will be a full matrix in general. However, for any factorization of the form B = SST$B=S{S}^{\mathrm{T}}$, the generalized symmetric problem reduces to the standard form
 S − 1AS − T(STx) = λ(STx) $S-1AS-T(STx)=λ(STx)$
and there does exist a factorization such that S1AST${S}^{-1}A{S}^{-\mathrm{T}}$ is still of band form (see Crawford (1973)). Writing
 C = S − 1AS − T  and  y = STx $C=S-1AS-T and y=STx$
the standard form is Cy = λy$Cy=\lambda y$ and the bandwidth of C$C$ is the maximum bandwidth of A$A$ and B$B$.
Each stage in the transformation consists of two phases. The first reduces a leading principal sub-matrix of B$B$ to the identity matrix and this introduces nonzero elements outside the band of A$A$. In the second, further transformations are applied which leave the reduced part of B$B$ unaltered and drive the extra elements upwards and off the top left corner of A$A$. Alternatively, B$B$ may be reduced to the identity matrix starting at the bottom right-hand corner and the extra elements introduced in A$A$ can be driven downwards.
The advantage of the ULDLTUT$ULD{L}^{\mathrm{T}}{U}^{\mathrm{T}}$ decomposition of B$B$ is that no extra elements have to be pushed over the whole length of A$A$. If k$k$ is taken as approximately n / 2$n/2$, the shifting is limited to halfway. At each stage the size of the triangular bumps produced in A$A$ depends on the number of rows and columns of B$B$ which are eliminated in the first phase and on the bandwidth of B$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$A$.
In this function, a is defined as being at least as wide as B$B$ and must be filled out with zeros if necessary as it is overwritten with C$C$. The number of rows and columns of B$B$ which are effectively eliminated at each stage is mA${m}_{A}$.

## Example

```function nag_matop_real_symm_posdef_geneig_example
k = int64(4);
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];
[b, ifail] = nag_matop_real_symm_posdef_fac(k, b);
[aOut, bOut, ifail] = nag_matop_real_symm_posdef_geneig(k, a, b)
```
```

aOut =

0    0.0685    0.1152    0.1325    0.1445    0.1563    0.1500    0.0952    0.0371
0.1692    0.0684    0.0207    0.0040   -0.0122   -0.0194    0.0479    0.1838    0.2898

bOut =

0    0.2178    0.2366    0.2460    0.2547    0.2636    0.2722    0.2792    0.2661
101.0000   97.2079   97.5581   91.7279   98.1475   98.6499   99.1822  100.2844  109.0000

ifail =

0

```
```function f01bv_example
k = int64(4);
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];
[b, ifail] = f01bu(k, b);
[aOut, bOut, ifail] = f01bv(k, a, b)
```
```

aOut =

0    0.0685    0.1152    0.1325    0.1445    0.1563    0.1500    0.0952    0.0371
0.1692    0.0684    0.0207    0.0040   -0.0122   -0.0194    0.0479    0.1838    0.2898

bOut =

0    0.2178    0.2366    0.2460    0.2547    0.2636    0.2722    0.2792    0.2661
101.0000   97.2079   97.5581   91.7279   98.1475   98.6499   99.1822  100.2844  109.0000

ifail =

0

```