- while (k <= max_iterations)
-
- if (norm(rk) < tolerance)
- % Check our stopping condition. This should catch the k=0 case.
- x = xk;
- return;
- end
-
+ while (k <= max_iterations && norm(rk, 'inf') > tolerance)
- alpha_k = rkzk/(dk' * Qdk);
+ % We're going to divide by this quantity...
+ dkQdk = dk' * Qdk;
+
+ % So if it's too close to zero, we replace it with something
+ % comparable but non-zero.
+ if (dkQdk < eps)
+ dkQdk = eps;
+ end
+
+ alpha_k = rkzk/dkQdk;
x_next = xk + (alpha_k * dk);
% The recursive definition of r_next is prone to accumulate
% roundoff error. When sqrt(n) divides k, we recompute the
x_next = xk + (alpha_k * dk);
% The recursive definition of r_next is prone to accumulate
% roundoff error. When sqrt(n) divides k, we recompute the