Replace step_length_cgm() with a direct call to step_length_positive_definite().

[octave.git] / optimization / vanilla_cgm.m

diff --git a/optimization/vanilla_cgm.m b/optimization/vanilla_cgm.m
index 2c94487a65508e876903343abf457cda853bb929..b511b484dfaab2ecba7fa6b8702be5080f977ba9 100644 (file)
--- a/optimization/vanilla_cgm.m
+++ b/optimization/vanilla_cgm.m
@@ -1,4 +1,4 @@
-function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
+function [x, k] = vanilla_cgm(A, b, x0, tolerance, max_iterations)
   %
   % Solve,
   %
@@ -36,28 +36,33 @@ function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
   %
   % All vectors are assumed to be *column* vectors.
   %
-  zero_vector = zeros(length(x0), 1);
+
+  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;
+  rk = A*xk - b; % The first residual must be computed the hard way.
   pk = -rk;
 
-  for k = [ 0 : max_iterations ]
-    if (norm(rk) < tolerance)
-       % Success.
-       return;
+  while (k <= max_iterations && norm(rk, 'inf') > tolerance)
+    alpha_k = step_length_positive_definite(rk, 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
 
-    alpha_k = step_length_cgm(rk, A, pk);
-    x_next = x + alpha_k*pk;
-    r_next = rk + alpha_k*A*pk;
     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