author Michael Orlitzky Fri, 22 Mar 2013 01:51:22 +0000 (21:51 -0400) committer Michael Orlitzky Fri, 22 Mar 2013 01:51:22 +0000 (21:51 -0400)

index a6401e5ec21ff963883368576deb59ff0e1b36d0..0acf5b27998012a85e733a872cb88411913de2af 100644 (file)
@@ -35,6 +35,9 @@ function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
%
% All vectors are assumed to be *column* vectors.
%
+  % The rather verbose name of this function was chosen to avoid
+  % conflicts with other implementations.
+  %
n = length(x0);
M = eye(n);

index 63943482c8c6dd1b47d916d9bdaf400521d6b992..eb2089f5b4e50d367c8237be9c166d5065f70982 100644 (file)
@@ -39,7 +39,7 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q,
%
% OUTPUT:
%
-  %   - ``x`` - The solution to Qx=b.
+  %   - ``x`` - The computed solution to Qx=b.
%
%   - ``k`` - The ending value of k; that is, the number of
%   iterations that were performed.
@@ -52,6 +52,9 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q,
% Conjugate-Gradient Method", we are supposed to define
% d_{0} = -z_{0}, not -r_{0} as written.
%
+  % The rather verbose name of this function was chosen to avoid
+  % conflicts with other implementations.
+  %
% REFERENCES:
%
%   1. Guler, Osman. Foundations of Optimization. New York, Springer,
@@ -59,7 +62,7 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q,
%

% Set k=0 first, that way the references to xk,rk,zk,dk which
-  % immediately follow correspond to x0,r0,z0,d0 respectively.
+  % immediately follow correspond (semantically) to x0,r0,z0,d0.
k = 0;

xk = x0;
@@ -68,9 +71,11 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q,
dk = -zk;

for k = [ 0 : max_iterations ]
+
if (norm(rk) < tolerance)
-       x = xk;
-       return;
+      % Check our stopping condition. This should catch the k=0 case.
+      x = xk;
+      return;
end

% Used twice, avoid recomputation.
@@ -80,6 +85,8 @@ function [x, k] = preconditioned_conjugate_gradient_method(Q,
% do them both, so we precompute the more expensive operation.
Qdk = Q * dk;

+    % After substituting the two previously-created variables, the
+    % following algorithm occurs verbatim in the reference.
alpha_k = rkzk/(dk' * Qdk);
x_next = xk + (alpha_k * dk);
r_next = rk + (alpha_k * Qdk);
@@ -94,5 +101,7 @@ 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.
x = xk;
end