function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
  %
  % Solve,
  %
  %   Ax = b
  %
  % or equivalently,
  %
  %   min [phi(x) = (1/2)*<Ax,x> + <b,x>]
  %
  % using the conjugate_gradient_method.
  %
  % INPUT:
  %
  %   - ``A`` -- The coefficient matrix of the system to solve. Must
  %     be positive definite.
  %
  %   - ``b`` -- The right-hand-side of the system to solve.
  %
  %   - ``x0`` -- The starting point for the search.
  %
  %   - ``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`` - 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.
  %
  % The rather verbose name of this function was chosen to avoid
  % conflicts with other implementations.
  %
  n = length(x0);
  M = eye(n);

  % 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
