function [x, k] = simple_preconditioned_cgm(Q, ...
					    M, ...
					    b, ...
					    x0, ...
					    tolerance, ...
					    max_iterations)
  %
  % Solve,
  %
  %   Qx = b
  %
  % or equivalently,
  %
  %   min [phi(x) = (1/2)*<Qx,x> + <b,x>]
  %
  % using the preconditioned conjugate gradient method (14.54 in
  % Guler).
  %
  % INPUT:
  %
  %   - ``Q`` -- The coefficient matrix of the system to solve. Must
  %     be positive definite.
  %
  %   - ``M`` -- The preconditioning matrix. If the actual matrix used
  %     to precondition ``Q`` is called ``C``, i.e. ``C^(-1) * Q *
  %     C^(-T) == \bar{Q}``, then M=CC^T. Must be symmetric positive-
  %     definite. See for example Golub and Van Loan.
  %
  %   - ``b`` -- The right-hand-side of the system to solve.
  %
  %   - ``x0`` -- The starting point for the search.
  %
  %   - ``tolerance`` -- How close ``Qx`` has to be to ``b`` (in
  %     magnitude) before we stop.
  %
  %   - ``max_iterations`` -- The maximum number of iterations to
  %     perform.
  %
  % OUTPUT:
  %
  %   - ``x`` - The solution to Qx=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.
  %
  % REFERENCES:
  %
  %   1. Guler, Osman. Foundations of Optimization. New York, Springer,
  %   2010.
  %

  % This isn't great in practice, since the CGM is usually used on
  % huge sparse systems.
  Ct = chol(M);
  C = Ct';
  C_inv = inv(C);
  Ct_inv = inv(Ct);

  Q_bar = C_inv * Q * Ct_inv;
  b_bar = C_inv * b;

  % But it sure is easy.
  [x_bar, k] = vanilla_cgm(Q_bar, b_bar, x0, tolerance, max_iterations);

  % The solution to Q_bar*x_bar == b_bar is x_bar = Ct*x.
  x = Ct_inv * x_bar;
end
