TQLI Procedure

Uses the QL algorithm with implicit shifts to determine the eigenvalues and eigenvectors of a real, symmetric, tridiagonal matrix.

TQLI, d, e, z

d An n-element vector. On input, it contains the diagonal elements of the matrix.

e An n-element vector. On input, it contains the off-diagonal elements of the matrix. e₀ is arbitrary.

d An n-element vector. On output, it contains the eigenvectors.

e An n-element vector. On output, it is destroyed.

z A matrix containing the n eigenvectors. The eigenvectors are stored by rows; for example, z(*, 0) contains the first eigenvector.

None.

The TRED2 procedure may be used to reduce a real symmetric matrix to the tridiagonal form that is suitable for input to TQLI.

If the eigenvectors of a tridiagonal matrix are desired, z should be input as an identity matrix. If the eigenvectors of a matrix that has been reduced by TRED2 are desired, z should be input as the matrix Q output by TRED2.

TQLI is based on the routine of the same name found in Numerical Recipes in C: The Art of Scientific Computing, by Flannery, Press, Teukolsky, and Vetterling, Cambridge University Press, Cambridge, MA, 1988. It is used by permission.

However, because the order of subscripts in PV-WAVE are reversed in comparison to those in Numerical Recipes, the z matrix is transposed from the result described in that book.

To determine the eigenvalues and eigenvectors of a real, symmetric matrix A:

asave = A

; Save original matrix, as TRED2 destroys its input.

TRED2, A, d, e

; Reduce matrix A to tridiagonal form.

TQLI, d, e, A

; Obtain n eigenvalues in d, and n eigenvectors in A(*, i).

To verify the operation of these routines, use the definition of eigenvalues and eigenvectors. Matrix A is said to have an eigenvector x and corresponding eigenvalue if:

A . x = x

This is demonstrated by the following code fragment:

FOR i=0, N_ELEMENTS(d)-1 DO BEGIN

; For each eigenvector/value:

PRINT, 'Eigenvalue', i, '=', d(i)

; Print the eigenvalue.

PRINT, 'Eigenvector = ', Aa(*,i)

; Print the computed eigenvector.

PRINT, 'A x / lambda = ', asave # A(*,i)/d(i)

; Print the eigenvector again. This row vector should be equal to the
; eigenvector printed above.

ENDFOR

TRED2