-function x_star = conjugate_gradient_method(A, b, x0, tolerance)
+function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
%
% Solve,
%
%
% min [phi(x) = (1/2)*<Ax,x> + <b,x>]
%
- % using Algorithm 5.2 in Nocedal and Wright.
+ % using the conjugate_gradient_method.
%
% INPUT:
%
% - ``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:
%
% All vectors are assumed to be *column* vectors.
%
- zero_vector = zeros(length(x0), 1);
-
- k = 0;
- xk = x0;
- rk = A*xk - b; % The first residual must be computed the hard way.
- pk = -rk;
-
- while (norm(rk) > tolerance)
- alpha_k = step_length_cgm(rk, A, pk);
- x_next = xk + 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;
- rk = r_next;
- pk = p_next;
- end
+ % The rather verbose name of this function was chosen to avoid
+ % conflicts with other implementations.
+ %
+ n = length(x0);
+ M = eye(n);
- x_star = xk;
+ % The standard CGM is equivalent to the preconditioned CGM is you
+ % use the identity matrix as your preconditioner.
+ [x, k] = preconditioned_conjugate_gradient_method(A,
+ M,
+ b,
+ x0,
+ tolerance,
+ max_iterations);
end