% ``x0`` -- An n-by-1 vector; an initial guess to the solution.
%
% ``tolerance`` -- (optional; default: 1e-10) the stopping tolerance.
- % we stop when the relative error (the norm of the
- % residual divided by the norm of ``b``) is less
- % than ``tolerance``.
+ % we stop when the relative error (the infinity norm
+ % of the residual divided by the infinity norm of
+ % ``b``) is less than ``tolerance``.
%
% ``max_iterations`` -- (optional; default: intmax) the maximum
% number of iterations we will perform.
%
% ``iterations`` -- The number of iterations taken.
%
- % ``residual_norms`` -- An n-by-iterations vector of the residual norms
- % at each iteration. If not requested, they will
- % not be computed to save space.
+ % ``residual_norms`` -- An n-by-iterations vector of the residual
+ % (infinity)norms at each iteration. If not
+ % requested, they will not be computed to
+ % save space.
%
- save_residual_norms = false;
- if (nargout > 2)
- save_residual_norms = true;
- residual_norms = [];
- end
if (nargin < 4)
tolerance = 1e-10;
max_iterations = intmax();
end
- % This creates a diagonal matrix from the diagonal entries of A.
- M = diag(diag(A));
- N = M - A;
-
- % Avoid recomputing this each iteration.
- b_norm = norm(b);
-
- k = 0;
- xk = x0;
- rk = b - A*xk;
-
- while (norm(rk) > tolerance*b_norm && k < max_iterations)
- x_next = M \ (N*xk + b);
- r_next = b - A*x_next;
-
- k = k + 1;
- xk = x_next;
- rk = r_next;
- if (save_residual_norms)
- residual_norms = [ residual_norms; norm(rk) ];
- end
+ M_of_A = @(A) diag(diag(A));
+ if (nargout > 2)
+ [x, iterations, residual_norms] = ...
+ classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
+ else
+ [x, iterations] = ...
+ classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
end
-
- iterations = k;
- x = xk;
end