function x = jacobi( A, b, xguess, tol, maxiter )
%x=jacobi( A, b, xguess, tol )
%
% Solves the linear system Ax=b iteratively to find x.
%
% INPUTS:
%   A      - The matrix of coefficients
%   b      - The rhs vector.
%   xguess - The initial guess for the solution.
%   tol    - The solution tolerance

% AUTHOR: James C. Sutherland
% Spring, 2008

% set up for iteration
r = abs(b-A*xguess);  % calculate the initial residual.
ad = diag(A);         % grab the diagonal of the A matrix.

niter=0;

while norm(r,inf)>tol
  x = xguess + 1./ad .* (b-A*xguess);
  xguess = x;
  r = b - A*x;
  niter = niter+1;
end

fprintf('Converged in %1.0f iterations.  Linf norm=%1.1e\n',niter,norm(r,inf));
