author Michael Orlitzky Fri, 22 Mar 2013 20:28:17 +0000 (16:28 -0400) committer Michael Orlitzky Fri, 22 Mar 2013 20:28:17 +0000 (16:28 -0400)
Make the roundoff error modification from the PCGM in the vanilla CGM.

@@ -36,26 +36,39 @@ function [x, k] = vanilla_cgm(A, b, x0, tolerance, max_iterations)
%
% All vectors are assumed to be *column* vectors.
%
+
+  sqrt_n = floor(sqrt(length(x0)));
+
k = 0;
-  x = x0; % Eschew the 'k' suffix on 'x' for simplicity.
-  rk = A*x - b; % The first residual must be computed the hard way.
+  xk = x0; % Eschew the 'k' suffix on 'x' for simplicity.
+  rk = A*xk - b; % The first residual must be computed the hard way.
pk = -rk;

-  for k = [ 0 : max_iterations ]
+  while (k <= max_iterations)
if (norm(rk) < tolerance)
% Success.
+       x = xk;
return;
end

alpha_k = step_length_cgm(rk, A, pk);
-    x_next = x + alpha_k*pk;
-    r_next = rk + alpha_k*A*pk;
+    x_next = xk + alpha_k*pk;
+
+    % Avoid accumulated roundoff errors.
+    if (mod(k, sqrt_n) == 0)
+      r_next = A*x_next - b;
+    else
+      r_next = rk + (alpha_k * A * pk);
+    end
+
beta_next = (r_next' * r_next)/(rk' * rk);
p_next = -r_next + beta_next*pk;

k = k + 1;
-    x = x_next;
+    xk = x_next;
rk = r_next;
pk = p_next;
end
+
+  x = xk;
end