author Michael Orlitzky Sat, 23 Mar 2013 02:30:26 +0000 (22:30 -0400) committer Michael Orlitzky Sat, 23 Mar 2013 02:30:26 +0000 (22:30 -0400)
Use the infinity norm (instead of the 2-norm) in the PCGM.
Simplify the stopping conditions in the PCGM.

index 35fecaa99cdb81b7f22836791797c0bb2fe1acac..cd886084a9766e1d785df59185c358ae69d75cc7 100644 (file)
@@ -77,14 +77,7 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q, ...
zk = M \ rk;
dk = -zk;

-  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)
% Used twice, avoid recomputation.
rkzk = rk' * zk;

@@ -92,13 +85,22 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q, ...
% do them both, so we precompute the more expensive operation.
Qdk = Q * dk;

-    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 (abs(dkQdk) < eps)
+      dkQdk = sign(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
-    % residual to minimize this error. This modification is due to the
-    % second reference.
+    % residual to minimize this error. This modification was suggested
+    % by the second reference.
if (mod(k, sqrt_n) == 0)
r_next = Q*x_next - b;
else
@@ -120,7 +122,6 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q, ...
dk = d_next;
end

-  % The algorithm didn't converge, but we still want to return the
-  % terminal value of xk.
+  % If we make it here, one of the two stopping conditions was met.
x = xk;
end