9b2ddac24377e5d02651f7b9e52a5ee3a4c0707c
1 function x_star = conjugate_gradient_method(A, b, x0, tolerance)
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 Algorithm 5.2 in Nocedal and Wright.
12 %
13 % INPUT:
14 %
15 % - ``A`` -- The coefficient matrix of the system to solve. Must
16 % be positive definite.
17 %
18 % - ``b`` -- The right-hand-side of the system to solve.
19 %
20 % - ``x0`` -- The starting point for the search.
21 %
22 % - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in
23 % magnitude) before we stop.
24 %
25 % OUTPUT:
26 %
27 % - ``x_star`` - The solution to Ax=b.
28 %
29 % NOTES:
30 %
31 % All vectors are assumed to be *column* vectors.
32 %
33 zero_vector = zeros(length(x0), 1);
35 k = 0;
36 xk = x0;
37 rk = A*xk - b; % The first residual must be computed the hard way.
38 pk = -rk;
40 while (norm(rk) > tolerance)
41 alpha_k = step_length_cgm(rk, A, pk);
42 x_next = xk + alpha_k*pk;
43 r_next = rk + alpha_k*A*pk;
44 beta_next = (r_next' * r_next)/(rk' * rk);
45 p_next = -r_next + beta_next*pk;
47 k = k + 1;
48 xk = x_next;
49 rk = r_next;
50 pk = p_next;
51 end
53 x_star = xk;
54 end