Page 139 - Compact Numerical Methods For Computers

128 Compact numerical methods for computers

algorithm 14 to perform the search are probably worthwhile. In a larger set of
timings run on both a Data General ECLIPSE and an IBM/370 model 168 in six
hexadecimal digit arithmetic, the ‘original Jacobi method was on average the
slowest of five methods tested. (Note that the convergence test at STEP 17 below can-
not be used.)

10.4. ORGANISATION OF THE JACOBI ALGORITHM
To reduce program code length, the following procedure performs left and right
multiplications on the matrix A (k) separately, rather than use the formulae
(10.33)-(10.37). This may cause some slight loss in accuracy of the computed
results (compared, for instance, to the procedure jacobi discussed in §10.5).
Convergence testing is a relatively complicated matter. In the following al-
gorithm, convergence is assessed by counting the number of rotations skipped
during one sweep of the matrix because the rotation angle is too small to have any
effect on the current matrix. Rotations needed to order the diagonal elements of
the matrix (hence the eigenvalues) are always performed. The test that the sine of
the rotation angle is smaller than the machine precision is used to decide (at STEP
10) if the rotation is to be carried out when the diagonal elements are in order.
Unfortunately, if two diagonal elements are close in value, then the rotation angle
may not be small even if the off-diagonal element to be set to zero is quite small, so
that it is unnecessary to perform the rotation. Thus at STEP 7, the algorithm tests to
see if the off-diagonal element has a magnitude which when multiplied by 100 is
incapable of adjusting the diagonal elements, in which case the rotation is skipped. A
safety check is provided on the number of sweeps which are carried out since it does
not seem possible to guarantee the convergence of the algorithm to the satisfaction of
the above tests, Even if the algorithm terminates after 30 sweeps (my arbitrary choice
for the safety check limit) the approximations to the eigenvalues and eigenvectors
may still be good, and it is recommended that they be checked by computing
residuals and other tests.

Algorithm 14. A Jacobi algorithm for eigensolutions of a real symmetric matrix
Procedure evJacobi(n: integer; {order of matrices}
var A,V : matrix; {matrix and eigenvector array}
hitev: boolean); {flag to initialize eigenvector
array to a unit matrix if TRUE}
{alg14.pas ==
a variant of the Jacobi algorithm for computing eigensolutions of a
real symmetric matrix
n is the order of the eigenproblem
A is the real symmetric matrix in full
V will be rotated by the Jacobi transformations.
initev is TRUE if we wish this procedure to initialize
V to a unit matrix before commencing; otherwise,
V is assumed to be initialized outside the procedure,
e.g. for computing the eigenvectors of a generalized
eigenproblem as in alg15.pas.
Copyright 1988 J.C.Nash