# NAG Library Routine Document

## 1Purpose

f02sdf finds the eigenvector corresponding to a given real eigenvalue for the generalized problem $Ax=\lambda Bx$, or for the standard problem $Ax=\lambda x$, where $A$ and $B$ are real band matrices.

## 2Specification

Fortran Interface
 Subroutine f02sdf ( n, ma1, mb1, a, lda, b, ldb, sym, rmu, vec, d, work,
 Integer, Intent (In) :: n, ma1, mb1, lda, ldb, lwork Integer, Intent (Inout) :: ifail Integer, Intent (Out) :: iwork(n) Real (Kind=nag_wp), Intent (In) :: relep, rmu Real (Kind=nag_wp), Intent (Inout) :: a(lda,n), b(ldb,n), d(30) Real (Kind=nag_wp), Intent (Out) :: vec(n), work(lwork) Logical, Intent (In) :: sym
#include nagmk26.h
 void f02sdf_ (const Integer *n, const Integer *ma1, const Integer *mb1, double a[], const Integer *lda, double b[], const Integer *ldb, const logical *sym, const double *relep, const double *rmu, double vec[], double d[], Integer iwork[], double work[], const Integer *lwork, Integer *ifail)

## 3Description

Given an approximation $\mu$ to a real eigenvalue $\lambda$ of the generalized eigenproblem $Ax=\lambda Bx$, f02sdf attempts to compute the corresponding eigenvector by inverse iteration.
f02sdf first computes lower and upper triangular factors, $L$ and $U$, of $A-\mu B$, using Gaussian elimination with interchanges, and then solves the equation $Ux=e$, where $e={\left(1,1,1,\dots ,1\right)}^{\mathrm{T}}$ – this is the first half iteration.
There are then three possible courses of action depending on the input value of ${\mathbf{d}}\left(1\right)$.
 1 ${\mathbf{d}}\left(1\right)=0$. This setting should be used if $\lambda$ is an ill-conditioned eigenvalue (provided the matrix elements do not vary widely in order of magnitude). In this case it is essential to accept only a vector found after one half iteration, and $\mu$ must be a very good approximation to $\lambda$. If acceptable growth is achieved in the solution of $Ux=e$, then the normalized $x$ is accepted as the eigenvector. If not, columns of an orthogonal matrix are tried in turn in place of $e$. If none of these give acceptable growth, the routine fails, indicating that $\mu$ was not a sufficiently good approximation to $\lambda$. 2 ${\mathbf{d}}\left(1\right)>0$. This setting should be used if $\mu$ is moderately close to an eigenvalue which is not ill-conditioned (provided the matrix elements do not differ widely in order of magnitude). If acceptable growth is achieved in the solution of $Ux=e$, the normalized $x$ is accepted as the eigenvector. If not, inverse iteration is performed. Up to $30$ iterations are allowed to achieve a vector and a correction to $\mu$ which together give acceptably small residuals. 3 ${\mathbf{d}}\left(1\right)<0$. This setting should be used if the elements of $A$ and $B$ vary widely in order of magnitude. Inverse iteration is performed, but a different convergence criterion is used.
See Section 9.3 for further details.
Note that the bandwidth of the matrix $A$ must not be less than the bandwidth of $B$. If this is not so, either $A$ must be filled out with zeros, or matrices $A$ and $B$ may be reversed and $1/\mu$ supplied as an approximation to the eigenvalue $1/\lambda$. Also it is assumed that $A$ and $B$ each have the same number of subdiagonals as superdiagonals. If this is not so, they must be filled out with zeros. If $A$ and $B$ are both symmetric, only the upper triangles need be supplied.
Peters G and Wilkinson J H (1979) Inverse iteration, ill-conditioned equations and Newton's method SIAM Rev. 21 339–360
Wilkinson J H (1965) The Algebraic Eigenvalue Problem Oxford University Press, Oxford
Wilkinson J H (1972) Inverse iteration in theory and practice Symposia Mathematica Volume X 361–379 Istituto Nazionale di Alta Matematica, Monograf, Bologna
Wilkinson J H (1974) Notes on inverse iteration and ill-conditioned eigensystems Acta Univ. Carolin. Math. Phys. 1–2 173–177
Wilkinson J H (1979) Kronecker's canonical form and the $QZ$ algorithm Linear Algebra Appl. 28 285–303

## 5Arguments

1:     $\mathbf{n}$ – IntegerInput
On entry: $n$, the order of the matrices $A$ and $B$.
Constraint: ${\mathbf{n}}\ge 1$.
2:     $\mathbf{ma1}$ – IntegerInput
On entry: the value ${m}_{A}+1$, where ${m}_{A}$ is the number of nonzero lines on each side of the diagonal of $A$. Thus the total bandwidth of $A$ is $2{m}_{A}+1$.
Constraint: $1\le {\mathbf{ma1}}\le {\mathbf{n}}$.
3:     $\mathbf{mb1}$ – IntegerInput
On entry: if ${\mathbf{mb1}}\le 0$, $B$ is assumed to be the unit matrix. Otherwise mb1 must specify the value ${m}_{B}+1$, where ${m}_{B}$ is the number of nonzero lines on each side of the diagonal of $B$. Thus the total bandwidth of $B$ is $2{m}_{B}+1$.
Constraint: ${\mathbf{mb1}}\le {\mathbf{ma1}}$.
4:     $\mathbf{a}\left({\mathbf{lda}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput/Output
On entry: the $n$ by $n$ band matrix $A$. The ${m}_{A}$ subdiagonals must be stored in the first ${m}_{A}$ rows of the array; the diagonal in the (${m}_{A}+1$)th row; and the ${m}_{A}$ superdiagonals in rows ${m}_{A}+2$ to $2{m}_{A}+1$. Each row of the matrix must be stored in the corresponding column of the array. For example, if $n=6$ and ${m}_{A}=2$ the storage scheme is:
 $* * a31 a42 a53 a64 * a21 a32 a43 a54 a65 a11 a22 a33 a44 a55 a66 a12 a23 a34 a45 a56 * a13 a24 a35 a46 * * .$
Elements of the array marked $*$ need not be set. The following code assigns the matrix elements within the band to the correct elements of the array:
```   Do 20 j = 1, n
Do 10 i = max(1,j-MA1+1), min(n,j+MA1-1)
a(i-j+MA1,j) = matrix(j,i)
10    Continue
20 Continue```
If ${\mathbf{sym}}=\mathrm{.TRUE.}$ (i.e., both $A$ and $B$ are symmetric), only the lower triangle of $A$ need be stored in the first ma1 rows of the array.
On exit: details of the factorization of $A-\stackrel{-}{\lambda }B$, where $\stackrel{-}{\lambda }$ is an estimate of the eigenvalue.
5:     $\mathbf{lda}$ – IntegerInput
On entry: the first dimension of the array a as declared in the (sub)program from which f02sdf is called.
Constraint: ${\mathbf{lda}}\ge 2×{\mathbf{ma1}}-1$.
6:     $\mathbf{b}\left({\mathbf{ldb}},{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput/Output
On entry: if ${\mathbf{mb1}}>0$, b must contain the $n$ by $n$ band matrix $B$, stored in the same way as $A$. If ${\mathbf{sym}}=\mathrm{.TRUE.}$, only the lower triangle of $B$ need be stored in the first mb1 rows of the array.
If ${\mathbf{mb1}}\le 0$, the array is not used.
On exit: elements in the top-left corner, and in the bottom right corner if ${\mathbf{sym}}=\mathrm{.FALSE.}$, are set to zero; otherwise the array is unchanged.
7:     $\mathbf{ldb}$ – IntegerInput
On entry: the first dimension of the array b as declared in the (sub)program from which f02sdf is called.
Constraints:
• if ${\mathbf{sym}}=\mathrm{.FALSE.}$, ${\mathbf{ldb}}\ge 2×{\mathbf{mb1}}-1$;
• if ${\mathbf{sym}}=\mathrm{.TRUE.}$, ${\mathbf{ldb}}\ge {\mathbf{mb1}}$.
8:     $\mathbf{sym}$ – LogicalInput
On entry: if ${\mathbf{sym}}=\mathrm{.TRUE.}$, both $A$ and $B$ are assumed to be symmetric and only their upper triangles need be stored. Otherwise sym must be set to .FALSE..
9:     $\mathbf{relep}$ – Real (Kind=nag_wp)Input
On entry: the relative error of the coefficients of the given matrices $A$ and $B$. If the value of relep is less than the machine precision, the machine precision is used instead.
10:   $\mathbf{rmu}$ – Real (Kind=nag_wp)Input
On entry: $\mu$, an approximation to the eigenvalue for which the corresponding eigenvector is required.
11:   $\mathbf{vec}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: the eigenvector, normalized so that the largest element is unity, corresponding to the improved eigenvalue ${\mathbf{rmu}}+{\mathbf{d}}\left(30\right)$.
12:   $\mathbf{d}\left(30\right)$ – Real (Kind=nag_wp) arrayInput/Output
On entry: ${\mathbf{d}}\left(1\right)$ must be set to indicate the type of problem (see Section 3):
${\mathbf{d}}\left(1\right)>0.0$
Indicates a well-conditioned eigenvalue.
${\mathbf{d}}\left(1\right)=0.0$
Indicates an ill-conditioned eigenvalue.
${\mathbf{d}}\left(1\right)<0.0$
Indicates that the matrices have elements varying widely in order of magnitude.
On exit: if ${\mathbf{d}}\left(1\right)\ne 0.0$ on entry, the successive corrections to $\mu$ are given in ${\mathbf{d}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,k$, where $k+1$ is the total number of iterations performed. The final correction is also given in the last position, ${\mathbf{d}}\left(30\right)$, of the array. The remaining elements of d are set to zero.
If ${\mathbf{d}}\left(1\right)=0.0$ on entry, no corrections to $\mu$ are computed and ${\mathbf{d}}\left(\mathit{i}\right)$ is set to $0.0$, for $\mathit{i}=1,2,\dots ,30$. Thus in all three cases the best available approximation to the eigenvalue is ${\mathbf{rmu}}+{\mathbf{d}}\left(30\right)$.
13:   $\mathbf{iwork}\left({\mathbf{n}}\right)$ – Integer arrayWorkspace
14:   $\mathbf{work}\left({\mathbf{lwork}}\right)$ – Real (Kind=nag_wp) arrayWorkspace
15:   $\mathbf{lwork}$ – IntegerInput
On entry: the dimension of the array work as declared in the (sub)program from which f02sdf is called.
Constraints:
• if ${\mathbf{d}}\left(1\right)\ne 0.0$, ${\mathbf{lwork}}\ge {\mathbf{n}}×\left({\mathbf{ma1}}+1\right)$;
• if ${\mathbf{d}}\left(1\right)=0.0$, ${\mathbf{lwork}}\ge 2×{\mathbf{n}}$.
16:   $\mathbf{ifail}$ – IntegerInput/Output
On entry: ifail must be set to $0$, $-1\text{​ or ​}1$. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{​ or ​}1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is $0$. When the value $-\mathbf{1}\text{​ or ​}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
 On entry, ${\mathbf{n}}<1$, or ${\mathbf{ma1}}<1$, or ${\mathbf{ma1}}>{\mathbf{n}}$, or ${\mathbf{lda}}<2×{\mathbf{ma1}}-1$, or ${\mathbf{ldb}}<{\mathbf{mb1}}$ when ${\mathbf{sym}}=\mathrm{.TRUE.}$, or ${\mathbf{ldb}}<2×{\mathbf{mb1}}-1$ when ${\mathbf{sym}}=\mathrm{.FALSE.}$ (ldb is not checked if ${\mathbf{mb1}}\le 0$).
${\mathbf{ifail}}=2$
 On entry, ${\mathbf{ma1}}<{\mathbf{mb1}}$. Either fill out a with zeros, or reverse the roles of a and b, and replace rmu by its reciprocal, i.e., solve $Bx={\lambda }^{-1}Ax\text{.}$
${\mathbf{ifail}}=3$
 On entry, ${\mathbf{lwork}}<2×{\mathbf{n}}$ when ${\mathbf{d}}\left(1\right)=0.0$, or ${\mathbf{lwork}}<{\mathbf{n}}×\left({\mathbf{ma1}}+1\right)$ when ${\mathbf{d}}\left(1\right)\ne 0.0$.
${\mathbf{ifail}}=4$
$A$ is null. If $B$ is nonsingular, all the eigenvalues are zero and any set of n orthogonal vectors forms the eigensolution.
${\mathbf{ifail}}=5$
$B$ is null. If $A$ is nonsingular, all the eigenvalues are infinite, and the columns of the unit matrix are eigenvectors.
${\mathbf{ifail}}=6$
 On entry, $A$ and $B$ are both null. The eigensolution is arbitrary.
${\mathbf{ifail}}=7$
${\mathbf{d}}\left(1\right)\ne 0.0$ on entry and convergence is not achieved in $30$ iterations. Either the eigenvalue is ill-conditioned or rmu is a poor approximation to the eigenvalue. See Section 9.3.
${\mathbf{ifail}}=8$
${\mathbf{d}}\left(1\right)=0.0$ on entry and no eigenvector has been found after $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{n}},5\right)$ back-substitutions. rmu is not a sufficiently good approximation to the eigenvalue.
${\mathbf{ifail}}=9$
${\mathbf{d}}\left(1\right)<0.0$ on entry and rmu is too inaccurate for the solution to converge.
${\mathbf{ifail}}=-99$
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

## 7Accuracy

The eigensolution is exact for some problem
 $A+Ex=μB+Fx,$
where $‖E‖,‖F‖$ are of the order of $\eta \left(‖A‖+\mu ‖B‖\right)$, where $\eta$ is the value used for relep.

## 8Parallelism and Performance

f02sdf is not threaded in any implementation.

### 9.1Timing

The time taken by f02sdf is approximately proportional to $n{\left(2{m}_{A}+1\right)}^{2}$ for factorization, and to $n\left(2{m}_{A}+1\right)$ for each iteration.

### 9.2Storage

The storage of the matrices $A$ and $B$ is designed for efficiency on a paged machine.
f02sdf will work with full matrices but it will do so inefficiently, particularly in respect of storage requirements.

### 9.3Algorithmic Details

Inverse iteration is performed according to the rule
 $A-μByr+1=Bxr$
 $xr+ 1=1αr+ 1yr+ 1$
where ${\alpha }_{r+1}$ is the element of ${y}_{r+1}$ of largest magnitude.
Thus:
 $A-μBxr+1=1αr+1Bxr.$
Hence the residual corresponding to ${x}_{r+1}$ is very small if $\left|{\alpha }_{r+1}\right|$ is very large (see Peters and Wilkinson (1979)). The first half iteration, $U{y}_{1}=e$, corresponds to taking ${L}^{-1}PB{x}_{0}=e$.
If $\mu$ is a very accurate eigenvalue, then there should always be an initial vector ${x}_{0}$ such that one half iteration gives a small residual and thus a good eigenvector. If the eigenvalue is ill-conditioned, then second and subsequent iterated vectors may not be even remotely close to an eigenvector of a neighbouring problem (see pages 374–376 of Wilkinson (1972) and Wilkinson (1974)). In this case it is essential to accept only a vector obtained after one half iteration.
However, for well-conditioned eigenvalues, there is no loss in performing more than one iteration (see page 376 of Wilkinson (1972)), and indeed it will be necessary to iterate if $\mu$ is not such a good approximation to the eigenvalue. When the iteration has converged, ${y}_{r+1}$ will be some multiple of ${x}_{r}$, ${y}_{r+1}={\beta }_{r+1}{x}_{r}$, say.
Therefore
 $A-μBβr+1xr=Bxr,$
giving
 $A-μ+1βr+ 1 B xr=0.$
Thus $\mu +\frac{1}{{\beta }_{r+1}}$ is a better approximation to the eigenvalue. ${\beta }_{r+1}$ is obtained as the element of ${y}_{r+1}$ which corresponds to the element of largest magnitude, $+1$, in ${x}_{r}$. The routine terminates when $‖\left(A-\left(\mu +\frac{1}{{\beta }_{r}}\right)B\right){x}_{r}‖$ is of the order of the machine precision relative to $‖A‖+\left|\mu \right|‖B‖$.
If the elements of $A$ and $B$ vary widely in order of magnitude, then $‖A‖$ and $‖B‖$ are excessively large and a different convergence test is required. The routine terminates when the difference between successive corrections to $\mu$ is small relative to $\mu$.
In practice one does not necessarily know if the given problem is well-conditioned or ill-conditioned. In order to provide some information on the condition of the eigenvalue or the accuracy of $\mu$ in the event of failure, successive values of $\frac{1}{{\beta }_{r}}$ are stored in the vector d when ${\mathbf{d}}\left(1\right)$ is nonzero on input. If these values appear to be converging steadily, then it is likely that $\mu$ was a poor approximation to the eigenvalue and it is worth trying again with ${\mathbf{rmu}}+{\mathbf{d}}\left(30\right)$ as the initial approximation. If the values in d vary considerably in magnitude, then the eigenvalue is ill-conditioned.
A discussion of the significance of the singularity of $A$ and/or $B$ is given in relation to the $QZ$ algorithm in Wilkinson (1979).

## 10Example

Given the generalized eigenproblem $Ax=\lambda Bx$ where
 $A= 1 1 2 -1 2 1 2 -1 3 1 2 -1 4 1 -1 5 and B= 5 1 1 4 2 2 3 2 2 2 1 1 1$
find the eigenvector corresponding to the approximate eigenvalue $-12.33$.
Although $B$ is symmetric, $A$ is not, so sym must be set to .FALSE. and all the elements of $B$ in the band must be supplied to the routine. $A$ (as written above) has $1$ subdiagonal and $2$ superdiagonals, so ma1 must be set to $3$ and $A$ filled out with an additional subdiagonal of zeros. Each row of the matrices is read in as data in turn.

### 10.1Program Text

Program Text (f02sdfe.f90)

### 10.2Program Data

Program Data (f02sdfe.d)

### 10.3Program Results

Program Results (f02sdfe.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017