1 function x_star = conjugate_gradient_method(A, b, x0, tolerance)
9 ## min [phi(x) = (1/2)*<Ax,x> + <b,x>]
11 ## using Algorithm 5.2 in Nocedal and Wright.
15 ## - ``A`` -- The coefficient matrix of the system to solve. Must
16 ## be positive definite.
18 ## - ``b`` -- The right-hand-side of the system to solve.
20 ## - ``x0`` -- The starting point for the search.
22 ## - ``tolerance`` -- How close ``Ax`` has to be to ``b`` (in
23 ## magnitude) before we stop.
27 ## - ``x_star`` - The solution to Ax=b.
31 ## All vectors are assumed to be *column* vectors.
33 zero_vector = zeros(length(x0), 1);
37 rk = A*xk - b; # The first residual must be computed the hard way.
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;