Page 148 - Compact Numerical Methods For Computers

The generalised symmetric matrix eigenvalue problem 137
plane rotations. These post-multiply the current approximation to the vectors,
which is usually initialised to the unit matrix l . From (11.11a), however, we have
n
- 1
X = ZD Y (11.20)
where Y is the matrix of eigenvectors of A , and hence a product of plane
2
rotations. Therefore, by setting the initial approximation of the vectors to ZD - 1
when solving the eigenproblem (11.10) of A , it is not necessary to solve (11.11a)
2
explicitly, nor to save Z or D.
The form (2.63) is not the only generalised eigenproblem which may arise.
Several are discussed in Wilkinson and Reinsch (1971, pp 303-14). This treat-
ment is based mainly on the Choleski decomposition. Nash (1974) also uses a
Choleski-type decomposition for the generalised eigenproblem (2.63) in the case
where A and B are complex Hermitian matrices. This can be solved by ‘doubling-
up’ the matrices to give the real symmetric eigenproblem

(11.21)

2
which therefore requires 12n matrix elements to take the expanded matrices and
resulting eigenvectors. Nash (1974) shows how it may be solved using only
2
4n + 4n matrix elements.

Algorithm 15. Solution of a generalised matrix eigenvalue problem by two applications
of the Jacobi algorithm

procedure genevJac( n : integer; {order of problem}
var A, B, V : rmatrix); {matrices and eigenvectors}
{alg15pas ==
To solve the generalized symmetric eigenvalue problem
Ax = eBx
where A and B are symmetric and B is positive-definite.
Method: Eigenvalue decomposition of B via Jacobi method to form the
‘matrix square root’ B-half. The Jacobi method is applied a second time
to the matrix
C = Bihalf A Bihalf
where Bihalf is the inverse of B-half. The eigenvectors x are the columns
of the matrix
X = Bihalf V
where V is the matrix whose columns are the eigenvectors of matrix C.
Copyright 1988 J.C.Nash
}
var
i,j,k,m : integer;
s : real;
initev : boolean;
begin {STEPS 0 and 1}
writeln(‘alg15.pas -- generalized symmetric matrix eigenproblem’);
initev := true; {to indicate eigenvectors must be initialized}
writeln(‘Eigensolutions of metric B’);
evJacobi(n, B, V, initev); {eigensolutions of B -- STEP 2}