Computes the eigenvalues of an n × n matrix, A. Optionally, the eigenvectors of A or AT are computed. Using the eigenvectors of A gives the decomposition AV = VE, where V is an n × n complex matrix of eigenvectors, and E is the complex diagonal matrix of eigenvalues. Other options include the reduction of A to upper triangular or Schur form, reduction to block upper triangular form with 2 × 2 or unit sized diagonal block matrices, and reduction to upper Hessenberg form.

The input matrix A is first balanced. The resulting similar matrix is transformed to upper Hessenberg form using orthogonal transformations. The double-shifted QR algorithm transforms the Hessenberg matrix so that 2 × 2 or unit sized blocks remain along the main diagonal. Any off-diagonal that is classified as “small” in order to achieve this block form is set to the value zero. Next the block upper triangular matrix is transformed to upper triangular form with unitary rotations. The eigenvectors of the upper triangular matrix are computed using back substitution. Care is taken to avoid overflows during this process. At the end, eigenvectors are normalized to have Euclidean length one, with the largest component real and positive. This algorithm follows that given in Golub and Van Loan, (1989, Chapter 7), with some novel organizational details for additional options, efficiency and robustness.

See the messages.gls file for error messages for LIN_EIG_GEN. These error messages are numbered 841–858; 861–878; 881–898; 901–918.

The eigenvalues of a random real matrix are computed. These values define a complex diagonal matrix E. Their correctness is checked by obtaining the eigenvector matrix V and verifying that the residuals R = AV ‑ VE are small. Also, see operator_ex29, supplied with the product examples.

The roots of a complex polynomial equation,

are required. This algebraic equation is formulated as a matrix eigenvalue problem. The equivalent matrix eigenvalue problem is solved using the upper Hessenberg matrix which has the value zero except in row number 1 and along the first subdiagonal. The entries in the first row are given by
a1,j = – bj, i = 1, …, n, while those on the first subdiagonal have the value one. This is a companion matrix for the polynomial. The results are checked by testing for small values of ∣f(ei)∣, i = 1, …, n, at the eigenvalues of the matrix, which are the roots of f(z). Also, see operator_ex30, supplied with the product examples.

The efficient solution of a family of linear algebraic equations is required. These systems are (A + hI)x = b. Here A is an n × n real matrix, I is the identity matrix, and b is the right-hand side matrix. The scalar h is such that the coefficient matrix is nonsingular. The method is based on the Schur form for matrix A: AW = WT, where W is unitary and T is upper triangular. This provides an efficient solution method for several values of h, once the Schur form is computed. The solution steps solve, for y, the upper triangular linear system

Then, x = x(h) = Wy. This is an efficient and accurate method for such parametric systems provided the expense of computing the Schur form has a pay-off in later efficiency. Using the Schur form in this way, it is not required to compute an LU factorization of A + hI with each new value of h. Note that even if the data A, h, and b are real, subexpressions for the solution may involve complex intermediate values, with x(h) finally a real quantity. Also, see operator_ex31, supplied with the product examples.

A matrix A has entries that are subject to uncertainty. This is expressed as the realization that A can be replaced by the matrix A + ηB, where the value η is “small” but still significantly larger than machine precision. The matrix B satisfies ∥B∥ ≤ ∥A∥. A variation in eigenvalues is estimated using analysis found in Golub and Van Loan, (1989, Chapter 7, p. 344). Each eigenvalue and eigenvector is expanded in a power series in η. With

and normalized eigenvectors, the bound

is satisfied. The vectors ui and vi are the ordinary and adjoint eigenvectors associated respectively with ei and its complex conjugate. This gives an upper bound on the size of the change to each ∣ei∣ due to changing the matrix data. The reciprocal

is defined as the condition number of ei. Also, see operator_ex32, Chapter 10.