iterative/classical_iteration.m

   1 function [x, iterations, residual_norms] = ...
   2          classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations)
   3
   4   save_residual_norms = false;
   5   if (nargout > 2)
   6      save_residual_norms = true;
   7      residual_norms = [];
   8   end
   9
  10   % This creates a diagonal matrix from the diagonal entries of A.
  11   M = M_of_A(A);
  12   N = M - A;
  13
  14   % Avoid recomputing this at the beginning of each iteration.
  15   b_norm = norm(b, 'inf');
  16   relative_tolerance = tolerance*b_norm;
  17
  18   k = 0;
  19   xk = x0;
  20   rk = b - A*xk;
  21   rk_norm = norm(rk, 'inf');
  22
  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;
  28
  29     k = k + 1;
  30     xk = x_next;
  31     rk = r_next;
  32
  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
  35     % iteration.
  36     rk_norm = norm(rk, 'inf');
  37     if (save_residual_norms)
  38       residual_norms = [ residual_norms; rk_norm ];
  39     end
  40   end
  41
  42    iterations = k;
  43    x = xk;
  44 end