Engineering and Scientific Subroutine Library for AIX Version 3 Release 3: Guide and Reference

SGEGV and DGEGV--Eigenvalues and, Optionally, the Eigenvectors of a Generalized Real Eigensystem, Az=wBz, where A and B Are Real General Matrices

These subroutines compute the eigenvalues and, optionally, the eigenvectors of a generalized real eigensystem, where A and B are real general matrices. Eigenvalues w are based on the two parts returned in vectors alpha and beta, such that w_i = alpha_i / beta_i for beta_i <> 0, and w_i = infinity for beta_i = 0. Eigenvectors are returned in matrix Z:

Az = wBz

Table 126. Data Types

Syntax

On Entry

iopt

indicates the type of computation to be performed, where:

If iopt = 0, eigenvalues only are computed.

If iopt = 1, eigenvalues and eigenvectors are computed.

Specified as: a fullword integer; iopt = 0 or 1.

a

is the real general matrix A of order n. Specified as: an lda by (at least) n array, containing numbers of the data type indicated in Table 126. On output, A is overwritten; that is, the original input is not preserved.

lda

is the leading dimension of the array specified for a. Specified as: a fullword integer; lda > 0 and lda >= n.

b

is the real general matrix B of order n. Specified as: an ldb by (at least) n array, containing numbers of the data type indicated in Table 126. On output, B is overwritten; that is, the original input is not preserved.

ldb

is the leading dimension of the array specified for b. Specified as: a fullword integer; ldb > 0 and ldb >= n.

alpha

See On Return.

beta

See On Return.

z

See On Return.

ldz

has the following meaning, where:

If iopt = 0, it is not used in the computation.

If iopt = 1, it is the leading dimension of the output array specified for z.

Specified as: a fullword integer. It must have the following value, where:

If iopt = 0, ldz > 0.

If iopt = 1, ldz > 0 and ldz >= n.

n

is the order of matrices A and B. Specified as: a fullword integer; n >= 0.

aux

has the following meaning:

If naux = 0 and error 2015 is unrecoverable, aux is ignored.

Otherwise, it is a storage work area used by this subroutine. Its size is specified by naux.

Specified as: an area of storage, containing numbers of the data type indicated in Table 126. On output, the contents are overwritten.

naux

is the size of the work area specified by aux--that is, the number of elements in aux. Specified as: a fullword integer, where:

If naux = 0 and error 2015 is unrecoverable, SGEGV and DGEGV dynamically allocate the work area used by the subroutine. The work area is deallocated before control is returned to the calling program.

Otherwise, naux >= 3n.

On Return

alpha

is the vector alpha of length n, containing the numerators of the eigenvalues of the generalized real eigensystem Az = wBz. Returned as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 126.

beta

is the vector beta of length n, containing the denominators of the eigenvalues of the generalized real eigensystem Az = wBz. Returned as: a one-dimensional array of (at least) length n, containing numbers of the data type indicated in Table 126.

z

has the following meaning, where:

If iopt = 0, it is not used in the computation.

If iopt = 1, it is the matrix Z of order n, containing the eigenvectors of the generalized real eigensystem, Az = wBz. The eigenvector in column i of matrix Z corresponds to the eigenvalue w_i, computed using the alpha_i and beta_i values. Each eigenvector is normalized so that the modulus of its largest element is 1.

Returned as: an ldz by (at least) n array, containing numbers of the data type indicated in Table 126.

Notes

1. When you specify iopt = 0, you must specify:
   • A positive value for ldz
   • A dummy argument for z (see Example 1)
2. Matrices A, B, and Z, vectors alpha and beta, and the work area specified for aux must have no common elements; otherwise, results are unpredictable. See Concepts.
3. You have the option of having the minimum required value for naux dynamically returned to your program. For details, see Using Auxiliary Storage in ESSL.

Function

The following steps describe the methods used to compute the eigenvalues and, optionally, the eigenvectors of a generalized real eigensystem, Az = wBz, where A and B are real general matrices. The methods are based upon Moler and Stewart's QZ algorithm. You must calculate the resulting eigenvalues w based on the two parts returned in vectors alpha and beta. Each eigenvalue is calculated as follows: w_i = alpha_i/beta_i for beta_i <> 0 and w_i = infinity for beta_i = 0. Eigenvalues are unordered, except that complex conjugate pairs appear consecutively with the eigenvalue having the positive imaginary part first.

• For iopt = 0, the eigenvalues are computed as follows:
  1. Simultaneously reduce A to upper Hessenberg form and B to upper triangular form using orthogonal transformations.
  2. Reduce A from upper Hessenberg form to quasi-upper triangular form while maintaining the upper triangular form of B using orthogonal transformations.
  3. Compute the eigenvalues of the generalized real eigensystem with A in quasi-upper triangular form and B in upper triangular form using orthogonal transformations.
  4. The numerators and denominators of the eigenvalues are returned in vectors alpha and beta, respectively.
• For iopt = 1, the eigenvalues and eigenvectors are computed as follows:
  1. Simultaneously reduce A to upper Hessenberg form and B to upper triangular form using and accumulating orthogonal transformations.
  2. Reduce A from upper Hessenberg form to quasi-upper triangular form while maintaining the upper triangular form of B using and accumulating orthogonal transformations.
  3. Compute the eigenvalues of the generalized real eigensystem with A in quasi-upper triangular form and B in upper triangular form using and accumulating orthogonal transformations.
  4. Compute the eigenvectors of the generalized real eigensystem with A in quasi-upper triangular form and B in upper triangular form using back substitution.
  5. The numerators and denominators of the eigenvalues are returned in vectors alpha and beta, respectively, and the eigenvectors are returned in matrix Z.

For more information on these methods, see references [39], [43], [58], [83], [63], [62], [87], [97], and [99]. If n is 0, no computation is performed. The results of the computations using short- and long-precision data can vary in accuracy.

These algorithms have a tendency to generate underflows that may hurt overall performance. The system default is to mask underflow, which improves the performance of these subroutines.

Error Conditions

Resource Errors

Error 2015 is unrecoverable, naux = 0, and unable to allocate work area.

Computational Errors

Eigenvalue (i) failed to converge after (xxx) iterations:

• The eigenvalues (w_j, j = i+1, i+2, ..., n) are correct.
• If iopt = 1, then Z is modified, but no eigenvectors are correct.
• A and B have been modified.
• The return code is set to 1.
• i and xxx can be determined at run time by use of the ESSL error-handling facilities. To obtain this information, you must use ERRSET to change the number of allowable errors for error code 2101 in the ESSL error option table; otherwise, the default value causes your program to terminate when this error occurs. See What Can You Do about ESSL Computational Errors?.

Input-Argument Errors

1. iopt <> 0 or 1
2. n < 0
3. lda <= 0
4. n > lda
5. ldb <= 0
6. n > ldb
7. ldz <= 0 and iopt = 1
8. n > ldz and iopt = 1
9. Error 2015 is recoverable or naux<>0, and naux is too small--that is, less than the minimum required value. Return code 2 is returned if error 2015 is recoverable.

Example 1

This example shows how to find the eigenvalues only of a real generalized eigensystem problem, AZ = wBZ, where:

• NAUX is equal to 3N.
• AUX contains 3N elements.
• LDZ is set to 1 to avoid an error condition.
• DUMMY is used as a placeholder for argument z, which is not used in the computation.
• On output, matrices A and B are overwritten.

Note:

These matrices are from page 257 in referenced text [62].

Call Statement and Input

           IOPT  A  LDA  B  LDB  ALPHA   BETA     Z    LDZ  N   AUX  NAUX
            |    |   |   |   |     |      |       |     |   |    |    |
CALL DGEGV( 0  , A , 3 , B , 3 , ALPHA , BETA , DUMMY , 1 , 3 , AUX , 9  )

        *                 *
        | 10.0  1.0   2.0 |
A    =  |  1.0  3.0  -1.0 |
        |  1.0  1.0   2.0 |
        *                 *

        *               *
        | 1.0  2.0  3.0 |
B    =  | 4.0  5.0  6.0 |
        | 7.0  8.0  9.0 |
        *               *

Output

         *                       *
         |  (4.778424, 0.000000) |
ALPHA =  | (-4.760580, 0.000000) |
         |  (2.769466, 0.000000) |
         *                       *

         *           *
         |  0.000000 |
BETA  =  |  0.934851 |
         | 15.446215 |
         *           *

Example 2

This example shows how to find the eigenvalues and eigenvectors of a real generalized eigensystem problem, AZ = wBZ, where:

• NAUX is equal to 3N.
• AUX contains 3N elements.
• On output, matrices A and B are overwritten.

Note:

These matrices are from page 263 in referenced text [62].

Call Statement and Input

           IOPT  A  LDA  B  LDB  ALPHA   BETA   Z  LDZ  N   AUX  NAUX
            |    |   |   |   |     |      |     |   |   |    |    |
CALL DGEGV( 1  , A , 5 , B , 5 , ALPHA , BETA , Z , 5 , 5 , AUX , 15 )

        *                          *
        | 2.0  3.0  4.0  5.0   6.0 |
        | 4.0  4.0  5.0  6.0   7.0 |
A    =  | 0.0  3.0  6.0  7.0   8.0 |
        | 0.0  0.0  2.0  8.0   9.0 |
        | 0.0  0.0  0.0  1.0  10.0 |
        *                          *

        *                             *
        | 1.0  -1.0  -1.0  -1.0  -1.0 |
        | 0.0   1.0  -1.0  -1.0  -1.0 |
B    =  | 0.0   0.0   1.0  -1.0  -1.0 |
        | 0.0   0.0   0.0   1.0  -1.0 |
        | 0.0   0.0   0.0   0.0   1.0 |
        *                             *

Output

         *                       *
         |  (7.950050, 0.000000) |
         | (-0.277338, 0.000000) |
ALPHA =  |  (2.149669, 0.000000) |
         |  (6.720718, 0.000000) |
         | (10.987556, 0.000000) |
         *                       *

         *          *
         | 0.374183 |
         | 1.480299 |
BETA  =  | 1.636872 |
         | 1.213574 |
         | 0.908837 |
         *          *

        *
        | (1.000000, 0.000000)  (-0.483408, 0.000000)   (0.540696, 0.000000)
        | (0.565497, 0.000000)   (1.000000, 0.000000)   (0.684441, 0.000000)
Z    =  | (0.180429, 0.000000)  (-0.661372, 0.000000)  (-1.000000, 0.000000)
        | (0.034182, 0.000000)   (0.180646, 0.000000)   (0.363671, 0.000000)
        | (0.003039, 0.000000)  (-0.017732, 0.000000)  (-0.041865, 0.000000)
        *
                                                                      *
                          (1.000000, 0.000000)  (-1.000000, 0.000000) |
                          (0.722065, 0.000000)  (-0.610415, 0.000000) |
                         (-0.089003, 0.000000)  (-0.116987, 0.000000) |
                         (-0.223599, 0.000000)   (0.038979, 0.000000) |
                          (0.050111, 0.000000)   (0.018653, 0.000000) |
                                                                      *