1 function [x, k] = vanilla_cgm(A, b, x0, tolerance, max_iterations)
2 %
3 % Solve,
4 %
5 % Ax = b
6 %
7 % or equivalently,
8 %
9 % min [phi(x) = (1/2)*<Ax,x> + <b,x>]
10 %
11 % using the conjugate_gradient_method (Algorithm 5.2 in Nocedal and
12 % Wright).
13 %
14 % INPUT:
15 %
16 % - ``A`` -- The coefficient matrix of the system to solve. Must
17 % be positive definite.
18 %
19 % - ``b`` -- The right-hand-side of the system to solve.
20 %
21 % - ``x0`` -- The starting point for the search.
22 %
23 % - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in
24 % magnitude) before we stop.
25 %
26 % - ``max_iterations`` -- The maximum number of iterations to perform.
27 %
28 % OUTPUT:
29 %
30 % - ``x`` - The solution to Ax=b.
31 %
32 % - ``k`` - The ending value of k; that is, the number of iterations that
33 % were performed.
34 %
35 % NOTES:
36 %
37 % All vectors are assumed to be *column* vectors.
38 %
40 sqrt_n = floor(sqrt(length(x0)));
42 k = 0;
43 xk = x0;
44 rk = A*xk - b; % The first residual must be computed the hard way.
45 pk = -rk;
47 while (k <= max_iterations && norm(rk, 'inf') > tolerance)
48 alpha_k = step_length_positive_definite(rk, A, pk);
49 x_next = xk + alpha_k*pk;
51 % Avoid accumulated roundoff errors.
52 if (mod(k, sqrt_n) == 0)
53 r_next = A*x_next - b;
54 else
55 r_next = rk + (alpha_k * A * pk);
56 end
58 beta_next = (r_next' * r_next)/(rk' * rk);
59 p_next = -r_next + beta_next*pk;
61 k = k + 1;
62 xk = x_next;
63 rk = r_next;
64 pk = p_next;
65 end
67 x = xk;
68 end