NAG Library Routine Document

f02sdf (withdraw_real_band_geneig)

1
Purpose

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

2
Specification

Fortran Interface
Subroutine f02sdf ( n, ma1, mb1, a, lda, b, ldb, sym, relep, rmu, vec, d, iwork, work, lwork, ifail)
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
C Header Interface
#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)

3
Description

Given an approximation μ to a real eigenvalue λ of the generalized eigenproblem Ax=λBx, f02sdf attempts to compute the corresponding eigenvector by inverse iteration.
f02sdf first computes lower and upper triangular factors, L and U, of A-μB, using Gaussian elimination with interchanges, and then solves the equation Ux=e, where e=1,1,1,,1T – this is the first half iteration.
There are then three possible courses of action depending on the input value of d1.
1. d1=0.
This setting should be used if λ 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 μ must be a very good approximation to λ. 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 μ was not a sufficiently good approximation to λ.
2. d1>0.
This setting should be used if μ 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 μ which together give acceptably small residuals.
3. d1<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/μ supplied as an approximation to the eigenvalue 1/λ. 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.

4
References

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

5
Arguments

1:     n – IntegerInput
On entry: n, the order of the matrices A and B.
Constraint: n1.
2:     ma1 – IntegerInput
On entry: the value mA+1, where mA is the number of nonzero lines on each side of the diagonal of A. Thus the total bandwidth of A is 2mA+1.
Constraint: 1ma1n.
3:     mb1 – IntegerInput
On entry: if mb10, B is assumed to be the unit matrix. Otherwise mb1 must specify the value mB+1, where mB is the number of nonzero lines on each side of the diagonal of B. Thus the total bandwidth of B is 2mB+1.
Constraint: mb1ma1.
4:     aldan – Real (Kind=nag_wp) arrayInput/Output
On entry: the n by n band matrix A. The mA subdiagonals must be stored in the first mA rows of the array; the diagonal in the (mA+1)th row; and the mA superdiagonals in rows mA+2 to 2mA+1. Each row of the matrix must be stored in the corresponding column of the array. For example, if n=6 and mA=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 sym=.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-λ-B, where λ- is an estimate of the eigenvalue.
5:     lda – IntegerInput
On entry: the first dimension of the array a as declared in the (sub)program from which f02sdf is called.
Constraint: lda2×ma1-1.
6:     bldbn – Real (Kind=nag_wp) arrayInput/Output
On entry: if mb1>0, b must contain the n by n band matrix B, stored in the same way as A. If sym=.TRUE., only the lower triangle of B need be stored in the first mb1 rows of the array.
If mb10, the array is not used.
On exit: elements in the top-left corner, and in the bottom right corner if sym=.FALSE., are set to zero; otherwise the array is unchanged.
7:     ldb – IntegerInput
On entry: the first dimension of the array b as declared in the (sub)program from which f02sdf is called.
Constraints:
  • if sym=.FALSE., ldb2×mb1-1;
  • if sym=.TRUE., ldbmb1.
8:     sym – LogicalInput
On entry: if sym=.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:     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:   rmu – Real (Kind=nag_wp)Input
On entry: μ, an approximation to the eigenvalue for which the corresponding eigenvector is required.
11:   vecn – Real (Kind=nag_wp) arrayOutput
On exit: the eigenvector, normalized so that the largest element is unity, corresponding to the improved eigenvalue rmu+d30.
12:   d30 – Real (Kind=nag_wp) arrayInput/Output
On entry: d1 must be set to indicate the type of problem (see Section 3):
d1>0.0
Indicates a well-conditioned eigenvalue.
d1=0.0
Indicates an ill-conditioned eigenvalue.
d1<0.0
Indicates that the matrices have elements varying widely in order of magnitude.
On exit: if d10.0 on entry, the successive corrections to μ are given in di, for i=1,2,,k, where k+1 is the total number of iterations performed. The final correction is also given in the last position, d30, of the array. The remaining elements of d are set to zero.
If d1=0.0 on entry, no corrections to μ are computed and di is set to 0.0, for i=1,2,,30. Thus in all three cases the best available approximation to the eigenvalue is rmu+d30.
13:   iworkn – Integer arrayWorkspace
14:   worklwork – Real (Kind=nag_wp) arrayWorkspace
15:   lwork – IntegerInput
On entry: the dimension of the array work as declared in the (sub)program from which f02sdf is called.
Constraints:
  • if d10.0, lworkn×ma1+1;
  • if d1=0.0, lwork2×n.
16:   ifail – IntegerInput/Output
On entry: ifail must be set to 0, -1 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 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 -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6
Error Indicators and Warnings

If on entry 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:
ifail=1
On entry,n<1,
orma1<1,
orma1>n,
orlda<2×ma1-1,
orldb<mb1 when sym=.TRUE.,
orldb<2×mb1-1 when sym=.FALSE. (ldb is not checked if mb10).
ifail=2
On entry,ma1<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=λ-1Ax.
ifail=3
On entry,lwork<2×n when d1=0.0,
orlwork<n×ma1+1 when d10.0.
ifail=4
A is null. If B is nonsingular, all the eigenvalues are zero and any set of n orthogonal vectors forms the eigensolution.
ifail=5
B is null. If A is nonsingular, all the eigenvalues are infinite, and the columns of the unit matrix are eigenvectors.
ifail=6
On entry,A and B are both null. The eigensolution is arbitrary.
ifail=7
d10.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.
ifail=8
d1=0.0 on entry and no eigenvector has been found after minn,5 back-substitutions. rmu is not a sufficiently good approximation to the eigenvalue.
ifail=9
d1<0.0 on entry and rmu is too inaccurate for the solution to converge.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
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.
ifail=-999
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7
Accuracy

The eigensolution is exact for some problem
A+Ex=μB+Fx,  
where E,F are of the order of ηA+μB, where η is the value used for relep.

8
Parallelism and Performance

f02sdf is not threaded in any implementation.

9
Further Comments

9.1
Timing

The time taken by f02sdf is approximately proportional to n 2mA+1 2 for factorization, and to n2mA+1 for each iteration.

9.2
Storage

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.3
Algorithmic Details

Inverse iteration is performed according to the rule
A-μByr+1=Bxr  
xr+ 1=1αr+ 1yr+ 1  
where αr+1 is the element of yr+1 of largest magnitude.
Thus:
A-μBxr+1=1αr+1Bxr.  
Hence the residual corresponding to xr+1 is very small if αr+1 is very large (see Peters and Wilkinson (1979)). The first half iteration, Uy1=e, corresponds to taking L-1PBx0=e.
If μ is a very accurate eigenvalue, then there should always be an initial vector x0 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 μ is not such a good approximation to the eigenvalue. When the iteration has converged, yr+1 will be some multiple of xr, yr+1=βr+1xr, say.
Therefore
A-μBβr+1xr=Bxr,  
giving
A-μ+1βr+ 1 B xr=0.  
Thus μ+ 1βr+1  is a better approximation to the eigenvalue. βr+1 is obtained as the element of yr+1 which corresponds to the element of largest magnitude, +1, in xr. The routine terminates when A-μ+ 1 βr B xr  is of the order of the machine precision relative to A+μ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 μ is small relative to μ.
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 μ in the event of failure, successive values of 1βr  are stored in the vector d when d1 is nonzero on input. If these values appear to be converging steadily, then it is likely that μ was a poor approximation to the eigenvalue and it is worth trying again with rmu+d30 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).

10
Example

Given the generalized eigenproblem Ax=λ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.1
Program Text

Program Text (f02sdfe.f90)

10.2
Program Data

Program Data (f02sdfe.d)

10.3
Program Results

Program Results (f02sdfe.r)