author Michael Orlitzky Tue, 19 Mar 2013 13:28:17 +0000 (09:28 -0400) committer Michael Orlitzky Tue, 19 Mar 2013 13:28:17 +0000 (09:28 -0400)

index 9b2ddac24377e5d02651f7b9e52a5ee3a4c0707c..2c94487a65508e876903343abf457cda853bb929 100644 (file)
@@ -1,4 +1,4 @@
-function x_star = conjugate_gradient_method(A, b, x0, tolerance)
+function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
%
% Solve,
%
@@ -8,7 +8,8 @@ function x_star = conjugate_gradient_method(A, b, x0, tolerance)
%
%   min [phi(x) = (1/2)*<Ax,x> + <b,x>]
%
-  % using Algorithm 5.2 in Nocedal and Wright.
+  % using the conjugate_gradient_method (Algorithm 5.2 in Nocedal and
+  % Wright).
%
% INPUT:
%
@@ -22,9 +23,14 @@ function x_star = conjugate_gradient_method(A, b, x0, tolerance)
%   - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in
%     magnitude) before we stop.
%
+  %   - ``max_iterations`` -- The maximum number of iterations to perform.
+  %
% OUTPUT:
%
-  %   - ``x_star`` - The solution to Ax=b.
+  %   - ``x`` - The solution to Ax=b.
+  %
+  %   - ``k`` - The ending value of k; that is, the number of iterations that
+  %     were performed.
%
% NOTES:
%
@@ -33,22 +39,25 @@ function x_star = conjugate_gradient_method(A, b, x0, tolerance)
zero_vector = zeros(length(x0), 1);

k = 0;
-  xk = x0;
-  rk = A*xk - b; % The first residual must be computed the hard way.
+  x = x0; % Eschew the 'k' suffix on 'x' for simplicity.
+  rk = A*x - b; % The first residual must be computed the hard way.
pk = -rk;

-  while (norm(rk) > tolerance)
+  for k = [ 0 : max_iterations ]
+    if (norm(rk) < tolerance)
+       % Success.
+       return;
+    end
+
alpha_k = step_length_cgm(rk, A, pk);
-    x_next = xk + alpha_k*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;
-    xk = x_next;
+    x = x_next;
rk = r_next;
pk = p_next;
end
-
-  x_star = xk;
end