Page 261 - Compact Numerical Methods For Computers

248 Compact numerical methods for computers
Algorithm 25. Rayleigh quotient minimisation by conjugate gradients (cont.)
to the size of the numbers in the problem.}
begin {conjugate gradients search for improved eigenvector}
{Now generate the first search direction.}
for i := 1 to n do t[i] := -g[i]; {STEP 8}
itn:= 0; {STEP 9}
repeat {Major cg loop}
itn := itn+l; {to count the conjugate gradient iterations}
matmul(n. A, t, yvec); {STEP 10}
matmul(n, B, t, zvec); ipr := ipr+l;
tat := 0.0; tbt := 0.0; xat := 0.0; xbt := 0.0; {STEP 11}
for i := 1 to n do
begin
xat := xat+X[i]*yvec[i]; tat := tat+t[i]*yvec[i];
xbt := xbt+X[i]*zvec[i]; tbt := tbt+t[i]*zvec[i];
end;
{STEP 12 -- formation and solution of quadratic equation}
u := tat*xbt-xat*tbt; v := tat*xbx-xax*tbt;
w := xat*xbx-xax*xbt; d := v*v-4.0*U*W; {the discriminant}
if d<0.0 then halt; {STEP 13 -- safety check}
{Note: we may want a more imaginative response to this result
of the computations. Here we will assume imaginary roots of
the quadradic cannot arise, but perform the check.}
d := sqrt(d); {Now compute the roots in a stable manner -- STEP 14}
if w0.0 then step := -2.0*w/(v+d) else step := 0.5*(d-v)/u;
{STEP 15 -- update vectors}
count := 0; {to count the number of unchanged vector components}
xax := 0.0; xbx := 0.0;
for i := l to n do
begin
avec[i] := avec[i]+step*yvec[i];
bvec[i] := bvec[i]+step*zvec[i];
w := X[i]; X[i] := w+step*t[i];
if (reltest+w)=(reltest+X[i]) then count := count+l;
xax := xax+X[i]*avec[i]; xbx := xbx+X[i]*bvec[i];
end; {loop on i}
if xbx<=tol then halt {to avoid zero divide if B singular}
else pn := xax/xbx; {STEP 16}
if (count<n) and (pn<rq) then {STEPS 17 & 18}
begin
rq := pn gg := 0.0; {STEP 19}
for i:= 1 to n do
begin
g[i] := 2.0*(avec[i]-pn*bvec[i])/xbx; gg := gg+g[i]*g[i];
end; {loop on i}
if gg>tol then {STEP 20}
begin {STEP 21}
xbt := 0.0; for i := 1 to n do xbt := xbt+X[i]*zvec[i];
{STEP 22 -- compute formula (19.52)}
tabt := 0.0; beta := 0.0;
for i := l to n do
begin
w := yvec[i]-pn*zvec[i]; tabt := tabt+t[i]*w;