Page 259 - Compact Numerical Methods For Computers

246 Compact numerical methods for computers
Substituting
q = - g (16.5)
and the Hessian (19.50) into the expression
T T
z=g Ht/t Ht (16.4)
for the parameter in the two-term conjugate gradient recurrence, and noting that
T
g t = 0 (19.51)
by virtue of the ‘exact’ linear searches, gives
T T T T
z=[g (A-RB)t- (g g)(x Bt)]/[t (A-RB) t]. (19.52)
The work of Geradin (1971) implies that z should be evaluated by computing the
inner products within the square brackets and subtracting. I prefer to perform the
subtraction within each element of the inner product to reduce the effect of digit
cancellation.
Finally, condition (19.51) permits the calculation of the new value of the
Rayleigh quotient to be simplified. Instead of the expression which results from
expanding (19.33), we have from (19.49) evaluted at (x+kt), with (19.51), the
expression
T T T T
R=(t Ax+kt At)/(t Bx+kt Bt). (19.53)
This expression is not used in the algorithm. Fried (1972) has suggested several
other formulae for the recurrence parameter z of equation (19.52). At the time of
writing, too little comparative testing has been carried out to suggest that one
such formula is superior to any other.
Algorithm 2.5. Rayleigh quotient minimisation by conjugate gradients

procedure rqmcg( n : integer; {order of matrices}
A, B : rmatrix; {matrices defining eigenproblem}
var X : rvector; {eigenvector approximation, on both
input and output to this procedure}
var ipr : integer; {on input, a limit to the number of
matrix products allowed, on output, the number of
matrix products used}
var rq : real); {Rayleigh quotient = eigenvalue approx.}
{alg25.pas == Rayleigh quotient minimization by conjugate gradients
Minimize Rayleigh quotient
X-transpose A X / X-transpose B X
thereby solving generalized symmetric matrix eigenproblem
A X = r q B X
for minimal eigenvalue rq and its associated eigenvector.
A and B are assumed symmetric, with B positive definite.
While we supply explicit matrices here, only matrix products
are needed of the form v = A u, w = B u.
Copyright 1988 J.C.Nash
}
var
count i, itn, itlimit : integer;