function [x,hist,r] = richardson(A,b,varargin)
% RICHARDSON The Richardson iterative method for a linear system
%
% x = richardson(A,b) returns the solution vector x to
% Ax=b computed via the Richardson method to a tolerance
% of 1e-8.
%
% [x,hist,r] = richardson(A,b,'key',value,...) defines optional
% arguments with the following keys
%
%   'tol' : the stopping tolerance, default 1e-8
%   'maxiter' : the maximum iteration, default 10000
%   'x0' : the initial solution, default x0 = zeros
%     If you provide a 
%   'verbose' : output intermediate data
%
% The output history is index by the number of matrix vector products
% performed in the iterative method.  If x0 is not the default, then there
% is an additional matrix-vector product formed to compute the initial
% residual.
%
% Example:
%   n = 100;
%   [A,xy] = randgeom(n);
%   L = diag(sum(A)) - A; % construct the Laplacian of A
%   S = speye(n) - 1/10*L;
%   b = zeros(n,1); b(1) = 1; b(end) = -1;
%   x = richardson(A,b);

% David F. Gleich
% Purdue University, Networks & Matrix Computations sample code
% 2011-09-07

% argument parsing
n = size(b,1);

args = struct('maxiter',10000,'tol',1e-8,'x0',[],'verbose',0);
if nargin>2
    argsu = struct(varargin{:}); % copy user supplied args
    for fn=fieldnames(args)'
        if isfield(argsu,fn{1}), args.(fn{1}) = argsu.(fn{1}); end
    end
end

if isempty(args.x0), x = zeros(n,1); r = b;
else x = args.x0; r = b-A*x; end

hist = zeros(args.maxiter,1);

for iter=1:args.maxiter
    x = x + r;
    r = b - A*x;
    hist(iter) = norm(r)/norm(b); % use relative residual
    if hist(iter) < args.tol, break; end % stop when done
    if args.verbose, fprintf('iter %7i relres %12e\n', iter, hist(iter)); end
end

hist = hist(1:iter); % only keep useful history

if hist(end) > args.tol, warning('richardson:didNotConverge',...
        'relative residual is %e after %i iterations', ...
        hist(end), args.maxiter); end