function [x, k] = vanilla_cgm(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 (Algorithm 5.2 in Nocedal and
  % Wright).
  %
  % 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.
  %

  sqrt_n = floor(sqrt(length(x0)));

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

  while (k <= max_iterations)
    if (norm(rk) < tolerance)
       % Success.
       x = xk;
       return;
    end

    alpha_k = step_length_cgm(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

    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

  x = xk;
end
