-function [x, k] = conjugate_gradient_method(A, b, x0, tolerance, max_iterations)
+function [x, k] = vanilla_cgm(A, b, x0, tolerance, max_iterations)
%
% Solve,
%
%
% All vectors are assumed to be *column* vectors.
%
- zero_vector = zeros(length(x0), 1);
+
+ sqrt_n = floor(sqrt(length(x0)));
k = 0;
- x = x0; % Eschew the 'k' suffix on 'x' for simplicity.
- rk = A*x - b; % The first residual must be computed the hard way.
+ xk = x0;
+ rk = A*xk - b; % The first residual must be computed the hard way.
pk = -rk;
- for k = [ 0 : max_iterations ]
- if (norm(rk) < tolerance)
- % Success.
- return;
+ while (k <= max_iterations && norm(rk, 'inf') > tolerance)
+ alpha_k = step_length_cgm(rk, A, pk);
+ x_next = xk + alpha_k*pk;
+
+ % Avoid accumulated roundoff errors.
+ if (mod(k, sqrt_n) == 0)
+ r_next = A*x_next - b;
+ else
+ r_next = rk + (alpha_k * A * pk);
end
- alpha_k = step_length_cgm(rk, A, pk);
- x_next = x + alpha_k*pk;
- r_next = rk + alpha_k*A*pk;
beta_next = (r_next' * r_next)/(rk' * rk);
p_next = -r_next + beta_next*pk;
k = k + 1;
- x = x_next;
+ xk = x_next;
rk = r_next;
pk = p_next;
end
+
+ x = xk;
end