1 function [x, iterations, residual_norms] = ...
2 classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations)
4 save_residual_norms = false;
6 save_residual_norms = true;
10 % This creates a diagonal matrix from the diagonal entries of A.
14 % Avoid recomputing this at the beginning of each iteration.
15 b_norm = norm(b, 'inf');
16 relative_tolerance = tolerance*b_norm;
21 rk_norm = norm(rk, 'inf');
23 while (rk_norm > relative_tolerance && k < max_iterations)
24 % This is almost certainly slower than updating the individual
25 % components of xk, but it's "fast enough."
26 x_next = M \ (N*xk + b);
27 r_next = b - A*x_next;
33 % We store the current residual norm so that, if we're saving
34 % them, we don't need to recompute it at the beginning of the next
36 rk_norm = norm(rk, 'inf');
37 if (save_residual_norms)
38 residual_norms = [ residual_norms; rk_norm ];