hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_eigen_real_band_geneig (f02sd)

Purpose

nag_eigen_real_band_geneig (f02sd) finds the eigenvector corresponding to a given real eigenvalue for the generalized problem Ax = λBxAx=λBx, or for the standard problem Ax = λxAx=λx, where AA and BB are real band matrices.

Syntax

[a, b, vec, d, ifail] = f02sd(ma1, mb1, a, b, sym, relep, rmu, d, 'n', n)
[a, b, vec, d, ifail] = nag_eigen_real_band_geneig(ma1, mb1, a, b, sym, relep, rmu, d, 'n', n)

Description

Given an approximation μμ to a real eigenvalue λλ of the generalized eigenproblem Ax = λBxAx=λBx, nag_eigen_real_band_geneig (f02sd) attempts to compute the corresponding eigenvector by inverse iteration.
nag_eigen_real_band_geneig (f02sd) first computes lower and upper triangular factors, LL and UU, of AμBA-μB, using Gaussian elimination with interchanges, and then solves the equation Ux = eUx=e, where e = (1,1,1,,1)Te=(1,1,1,,1)T – this is the first half iteration.
There are then three possible courses of action depending on the input value of d(1)d1.
  1. d(1) = 0d1=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 = eUx=e, then the normalized xx is accepted as the eigenvector. If not, columns of an orthogonal matrix are tried in turn in place of ee. If none of these give acceptable growth, the function fails, indicating that μμ was not a sufficiently good approximation to λλ.
  2. d(1) > 0d1>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 = eUx=e, the normalized xx is accepted as the eigenvector. If not, inverse iteration is performed. Up to 3030 iterations are allowed to achieve a vector and a correction to μμ which together give acceptably small residuals.
  3. d(1) < 0d1<0.
    This setting should be used if the elements of AA and BB vary widely in order of magnitude. Inverse iteration is performed, but a different convergence criterion is used.
See Section [Algorithmic Details] for further details.
Note that the bandwidth of the matrix AA must not be less than the bandwidth of BB. If this is not so, either AA must be filled out with zeros, or matrices AA and BB may be reversed and 1 / μ1/μ supplied as an approximation to the eigenvalue 1 / λ1/λ. Also it is assumed that AA and BB each have the same number of subdiagonals as superdiagonals. If this is not so, they must be filled out with zeros. If AA and BB are both symmetric, only the upper triangles need be supplied.

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 QZQZ algorithm Linear Algebra Appl. 28 285–303

Parameters

Compulsory Input Parameters

1:     ma1 – int64int32nag_int scalar
The value mA + 1mA+1, where mAmA is the number of nonzero lines on each side of the diagonal of AA. Thus the total bandwidth of AA is 2mA + 12mA+1.
Constraint: 1ma1n1ma1n.
2:     mb1 – int64int32nag_int scalar
If mb10mb10, BB is assumed to be the unit matrix. Otherwise mb1 must specify the value mB + 1mB+1, where mBmB is the number of nonzero lines on each side of the diagonal of BB. Thus the total bandwidth of BB is 2mB + 12mB+1.
Constraint: mb1ma1mb1ma1.
3:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint lda2 × ma11lda2×ma1-1.
The nn by nn band matrix AA. The mAmA subdiagonals must be stored in the first mAmA rows of the array; the diagonal in the (mA + 1mA+1)th row; and the mAmA superdiagonals in rows mA + 2mA+2 to 2mA + 12mA+1. Each row of the matrix must be stored in the corresponding column of the array. For example, if n = 6n=6 and mA = 2mA=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 * *
.
* * 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:
 for j=1:n for i=max(1,j-ma1+1):min(n,j+ma1-1) a(i-j+ma1,j) = matrix(i,j); end end 
If sym = truesym=true (i.e., both AA and BB are symmetric), only the lower triangle of AA need be stored in the first ma1 rows of the array.
4:     b(ldb,n) – double array
ldb, the first dimension of the array, must satisfy the constraint
  • if sym = falsesym=false, ldb2 × mb11ldb2×mb1-1;
  • if sym = truesym=true, ldbmb1ldbmb1.
If mb1 > 0mb1>0, b must contain the nn by nn band matrix BB, stored in the same way as AA. If sym = truesym=true, only the lower triangle of BB need be stored in the first mb1 rows of the array.
If mb10mb10, the array is not used.
5:     sym – logical scalar
If sym = truesym=true, both AA and BB are assumed to be symmetric and only their upper triangles need be stored. Otherwise sym must be set to false.
6:     relep – double scalar
The relative error of the coefficients of the given matrices AA and BB. If the value of relep is less than the machine precision, the machine precision is used instead.
7:     rmu – double scalar
μμ, an approximation to the eigenvalue for which the corresponding eigenvector is required.
8:     d(3030) – double array
d(1)d1 must be set to indicate the type of problem (see Section [Description]):
d(1) > 0.0d1>0.0
Indicates a well-conditioned eigenvalue.
d(1) = 0.0d1=0.0
Indicates an ill-conditioned eigenvalue.
d(1) < 0.0d1<0.0
Indicates that the matrices have elements varying widely in order of magnitude.

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.)
nn, the order of the matrices AA and BB.
Constraint: n1n1.

Input Parameters Omitted from the MATLAB Interface

lda ldb iwork work lwork

Output Parameters

1:     a(lda,n) – double array
lda2 × ma11lda2×ma1-1.
Details of the factorization of AλBA-λ-B, where λλ- is an estimate of the eigenvalue.
2:     b(ldb,n) – double array
Elements in the top-left corner, and in the bottom right corner if sym = falsesym=false, are set to zero; otherwise the array is unchanged.
3:     vec(n) – double array
The eigenvector, normalized so that the largest element is unity, corresponding to the improved eigenvalue rmu + d(30)rmu+d30.
4:     d(3030) – double array
If d(1)0.0d10.0 on entry, the successive corrections to μμ are given in d(i)di, for i = 1,2,,ki=1,2,,k, where k + 1k+1 is the total number of iterations performed. The final correction is also given in the last position, d(30)d30, of the array. The remaining elements of d are set to zero.
If d(1) = 0.0d1=0.0 on entry, no corrections to μμ are computed and d(i)di is set to 0.00.0, for i = 1,2,,30i=1,2,,30. Thus in all three cases the best available approximation to the eigenvalue is rmu + d(30)rmu+d30.
5:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,n < 1n<1,
orma1 < 1ma1<1,
orma1 > nma1>n,
orlda < 2 × ma11lda<2×ma1-1,
orldb < mb1ldb<mb1 when sym = truesym=true,
orldb < 2 × mb11ldb<2×mb1-1 when sym = falsesym=false (ldb is not checked if mb10mb10).
  ifail = 2ifail=2
On entry,ma1 < mb1ma1<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.Bx=λ-1Ax.
  ifail = 3ifail=3
On entry,lwork < 2 × nlwork<2×n when d(1) = 0.0d1=0.0,
orlwork < n × (ma1 + 1)lwork<n×(ma1+1) when d(1)0.0d10.0.
  ifail = 4ifail=4
AA is null. If BB is nonsingular, all the eigenvalues are zero and any set of n orthogonal vectors forms the eigensolution.
  ifail = 5ifail=5
BB is null. If AA is nonsingular, all the eigenvalues are infinite, and the columns of the unit matrix are eigenvectors.
  ifail = 6ifail=6
On entry,AA and BB are both null. The eigensolution is arbitrary.
  ifail = 7ifail=7
d(1)0.0d10.0 on entry and convergence is not achieved in 3030 iterations. Either the eigenvalue is ill-conditioned or rmu is a poor approximation to the eigenvalue. See Section [Algorithmic Details].
  ifail = 8ifail=8
d(1) = 0.0d1=0.0 on entry and no eigenvector has been found after min (n,5)min(n,5) back-substitutions. rmu is not a sufficiently good approximation to the eigenvalue.
  ifail = 9ifail=9
d(1) < 0.0d1<0.0 on entry and rmu is too inaccurate for the solution to converge.

Accuracy

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

Further Comments

Timing

The time taken by nag_eigen_real_band_geneig (f02sd) is approximately proportional to n(2mA + 1)2n (2mA+1) 2 for factorization, and to n(2mA + 1)n(2mA+1) for each iteration.

Storage

The storage of the matrices AA and BB is designed for efficiency on a paged machine.
nag_eigen_real_band_geneig (f02sd) will work with full matrices but it will do so inefficiently, particularly in respect of storage requirements.

Algorithmic Details

Inverse iteration is performed according to the rule
(AμB)yr + 1 = Bxr
(A-μB)yr+1=Bxr
xr + 1 = 1/(αr + 1)yr + 1
xr+ 1=1αr+ 1yr+ 1
where αr + 1αr+1 is the element of yr + 1yr+1 of largest magnitude.
Thus:
(AμB)xr + 1 = 1/(αr + 1)Bxr.
(A-μB)xr+1=1αr+1Bxr.
Hence the residual corresponding to xr + 1xr+1 is very small if |αr + 1||αr+1| is very large (see Peters and Wilkinson (1979)). The first half iteration, Uy1 = eUy1=e, corresponds to taking L1PBx0 = eL-1PBx0=e.
If μμ is a very accurate eigenvalue, then there should always be an initial vector x0x0 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 + 1yr+1 will be some multiple of xrxr, yr + 1 = βr + 1xryr+1=βr+1xr, say.
Therefore
(AμB)βr + 1xr = Bxr,
(A-μB)βr+1xr=Bxr,
giving
(A(μ + 1/(βr + 1))B) xr = 0.
(A-(μ+1βr+ 1 ) B) xr=0.
Thus μ + 1/(βr + 1) μ+ 1βr+1  is a better approximation to the eigenvalue. βr + 1βr+1 is obtained as the element of yr + 1yr+1 which corresponds to the element of largest magnitude, + 1+1, in xrxr. The function terminates when (A(μ + 1/(βr))B)xr (A-(μ+ 1 βr ) B ) xr  is of the order of the machine precision relative to A + |μ|BA+|μ|B.
If the elements of AA and BB vary widely in order of magnitude, then AA and BB are excessively large and a different convergence test is required. The function 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) 1βr  are stored in the vector d when d(1)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 + d(30)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 AA and/or BB is given in relation to the QZQZ algorithm in Wilkinson (1979).

Example

function nag_eigen_real_band_geneig_example
ma1 = int64(3);
mb1 = int64(2);
a = [0, 0, 0, 0, 0;
     0, -1, -1, -1, -1;
     1, 2, 3, 4, 5;
     1, 1, 1, 1, 0;
     2, 2, 2, 0, 0];
b = [0, 1, 2, 2, 1;
     5, 4, 3, 2, 1;
     1, 2, 2, 1, 0];
sym = false;
relep = 0;
rmu = -12.33;
d = zeros(30, 1);
d(1) = 1;
[aOut, bOut, vec, dOut, ifail] = nag_eigen_real_band_geneig(ma1, mb1, a, b, sym, relep, rmu, d)
 

aOut =

    0.0160    0.0204    0.0423    0.0883   68.1366
   13.3300   25.2983   28.6600   17.3300         0
    2.0000    2.0000   13.3300         0         0
         0         0         0         0         0
         0         0         0         0         0


bOut =

     0     1     2     2     1
     5     4     3     2     1
     1     2     2     1     0


vec =

   -0.0572
    0.3951
   -0.8427
    1.0000
   -0.6540


dOut =

   -0.0094
   -0.0094
   -0.0094
   -0.0094
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
   -0.0094


ifail =

                    0


function f02sd_example
ma1 = int64(3);
mb1 = int64(2);
a = [0, 0, 0, 0, 0;
     0, -1, -1, -1, -1;
     1, 2, 3, 4, 5;
     1, 1, 1, 1, 0;
     2, 2, 2, 0, 0];
b = [0, 1, 2, 2, 1;
     5, 4, 3, 2, 1;
     1, 2, 2, 1, 0];
sym = false;
relep = 0;
rmu = -12.33;
d = zeros(30, 1);
d(1) = 1;
[aOut, bOut, vec, dOut, ifail] = f02sd(ma1, mb1, a, b, sym, relep, rmu, d)
 

aOut =

    0.0160    0.0204    0.0423    0.0883   68.1366
   13.3300   25.2983   28.6600   17.3300         0
    2.0000    2.0000   13.3300         0         0
         0         0         0         0         0
         0         0         0         0         0


bOut =

     0     1     2     2     1
     5     4     3     2     1
     1     2     2     1     0


vec =

   -0.0572
    0.3951
   -0.8427
    1.0000
   -0.6540


dOut =

   -0.0094
   -0.0094
   -0.0094
   -0.0094
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
   -0.0094


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013