From 8df0aa3b8e47b596626cf0f833d9e0c0143bf4d5 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Fri, 22 Mar 2013 16:28:17 -0400 Subject: [PATCH] Clean up the loop in the vanilla CGM. Make the roundoff error modification from the PCGM in the vanilla CGM. --- optimization/vanilla_cgm.m | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/optimization/vanilla_cgm.m b/optimization/vanilla_cgm.m index 3a72b14..9ada2cb 100644 --- a/optimization/vanilla_cgm.m +++ b/optimization/vanilla_cgm.m @@ -36,26 +36,39 @@ function [x, k] = vanilla_cgm(A, b, x0, tolerance, max_iterations) % % All vectors are assumed to be *column* vectors. % + + 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; % Eschew the 'k' suffix on 'x' for simplicity. + rk = A*xk - b; % The first residual must be computed the hard way. pk = -rk; - for k = [ 0 : max_iterations ] + while (k <= max_iterations) if (norm(rk) < tolerance) % Success. + x = xk; return; end alpha_k = step_length_cgm(rk, A, pk); - x_next = x + alpha_k*pk; - r_next = rk + alpha_k*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 + 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 -- 2.44.2