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_matop_real_symm_posdef_geneig (f01bv)

Purpose

nag_matop_real_symm_posdef_geneig (f01bv) transforms the generalized symmetric-definite eigenproblem Ax = λbxAx=λbx to the equivalent standard eigenproblem Cy = λyCy=λy, where AA, bb and CC are symmetric band matrices and bb is positive definite. bb 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

AA is a symmetric band matrix of order nn and bandwidth 2mA + 12mA+1. The positive definite symmetric band matrix BB, of order nn and bandwidth 2mB + 12mB+1, must have been previously decomposed by nag_matop_real_symm_posdef_fac (f01bu) as ULDLTUTULDLTUT. nag_matop_real_symm_posdef_geneig (f01bv) applies UU, LL and DD to AA, mAmA rows at a time, restoring the band form of AA at each stage by plane rotations. The parameter kk defines the change-over point in the decomposition of BB 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, kk should be chosen to be the multiple of mAmA nearest to n / 2n/2. The resulting symmetric band matrix CC is overwritten on a. The eigenvalues of CC, 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
kk, 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 BB.
Constraint: mb11knmb1-1kn.
2:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldama1ldama1.
The upper triangle of the nn by nn symmetric band matrix AA, with the diagonal of the matrix stored in the (mA + 1)(mA+1)th row of the array, and the mAmA superdiagonals within the band stored in the first mAmA rows of the array. Each column of the matrix is stored in the corresponding column of the array. For example, if n = 6n=6 and mA = 2mA=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 ldbmb1ldbmb1.
The elements of the decomposition of matrix BB 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.)
nn, the order of the matrices AA, BB and CC.
2:     ma1 – int64int32nag_int scalar
Default: The first dimension of the array a.
mA + 1mA+1, where mAmA is the number of nonzero superdiagonals in AA. Normally ma1nma1n.
3:     mb1 – int64int32nag_int scalar
Default: The first dimension of the array b.
mB + 1mB+1, where mBmB is the number of nonzero superdiagonals in BB.
Constraint: mb1ma1mb1ma1.

Input Parameters Omitted from the MATLAB Interface

m3 lda ldb v ldv w

Output Parameters

1:     a(lda,n) – double array
ldama1ldama1.
Stores the corresponding elements of CC.
2:     b(ldb,n) – double array
ldbmb1ldbmb1.
The elements of b will have been permuted.
3:     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,mb1 > ma1mb1>ma1.

Accuracy

In general the computed system is exactly congruent to a problem (A + E)x = λ(B + F)x(A+E)x=λ(B+F)x, where EE and FF are of the order of εκ(B)Aεκ(B)A and εκ(B)Bεκ(B)B respectively, where κ(B)κ(B) is the condition number of BB with respect to inversion and εε is the machine precision. This means that when BB is positive definite but not well-conditioned with respect to inversion, the method, which effectively involves the inversion of BB, may lead to a severe loss of accuracy in well-conditioned eigenvalues.

Further Comments

The time taken by nag_matop_real_symm_posdef_geneig (f01bv) is approximately proportional to n2mB2n2mB2 and the distance of kk from n / 2n/2, e.g., k = n / 4k=n/4 and k = 3n / 4k=3n/4 take 502%502% longer.
When BB is positive definite and well-conditioned with respect to inversion, the generalized symmetric eigenproblem can be reduced to the standard symmetric problem Py = λyPy=λy where P = L1ALTP=L-1AL-T and B = LLTB=LLT, the Cholesky factorization.
When AA and BB 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 PP since it will be a full matrix in general. However, for any factorization of the form B = SSTB=SST, the generalized symmetric problem reduces to the standard form
S1AST(STx) = λ(STx)
S-1AS-T(STx)=λ(STx)
and there does exist a factorization such that S1ASTS-1AS-T is still of band form (see Crawford (1973)). Writing
C = S1AST  and  y = STx
C=S-1AS-T  and  y=STx
the standard form is Cy = λyCy=λy and the bandwidth of CC is the maximum bandwidth of AA and BB.
Each stage in the transformation consists of two phases. The first reduces a leading principal sub-matrix of BB to the identity matrix and this introduces nonzero elements outside the band of AA. In the second, further transformations are applied which leave the reduced part of BB unaltered and drive the extra elements upwards and off the top left corner of AA. Alternatively, BB may be reduced to the identity matrix starting at the bottom right-hand corner and the extra elements introduced in AA can be driven downwards.
The advantage of the ULDLTUTULDLTUT decomposition of BB is that no extra elements have to be pushed over the whole length of AA. If kk is taken as approximately n / 2n/2, the shifting is limited to halfway. At each stage the size of the triangular bumps produced in AA depends on the number of rows and columns of BB which are eliminated in the first phase and on the bandwidth of BB. 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 AA.
In this function, a is defined as being at least as wide as BB and must be filled out with zeros if necessary as it is overwritten with CC. The number of rows and columns of BB which are effectively eliminated at each stage is mAmA.

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



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