1 function [x, k] = vanilla_cgm(A, b, x0, tolerance, max_iterations)
9 % min [phi(x) = (1/2)*<Ax,x> + <b,x>]
11 % using the conjugate_gradient_method (Algorithm 5.2 in Nocedal and
16 % - ``A`` -- The coefficient matrix of the system to solve. Must
17 % be positive definite.
19 % - ``b`` -- The right-hand-side of the system to solve.
21 % - ``x0`` -- The starting point for the search.
23 % - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in
24 % magnitude) before we stop.
26 % - ``max_iterations`` -- The maximum number of iterations to perform.
30 % - ``x`` - The solution to Ax=b.
32 % - ``k`` - The ending value of k; that is, the number of iterations that
37 % All vectors are assumed to be *column* vectors.
40 sqrt_n = floor(sqrt(length(x0)));
43 xk = x0; % Eschew the 'k' suffix on 'x' for simplicity.
44 rk = A*xk - b; % The first residual must be computed the hard way.
47 while (k <= max_iterations)
48 if (norm(rk) < tolerance)
54 alpha_k = step_length_cgm(rk, A, pk);
55 x_next = xk + alpha_k*pk;
57 % Avoid accumulated roundoff errors.
58 if (mod(k, sqrt_n) == 0)
59 r_next = A*x_next - b;
61 r_next = rk + (alpha_k * A * pk);
64 beta_next = (r_next' * r_next)/(rk' * rk);
65 p_next = -r_next + beta_next*pk;