function x=jacobi(A,b,x,M,tau);
% JACOBI: Solves Ax=b using Jacobi iterations
% USAGE:
%   x=JACOBI(A,b,x,M,tau,lambda);
%
% INPUTS:
%   A : nxn matrix
%   b : nxm matrix
%  [x] : nxm matrix of starting values. Default: b
%  [M] : maximum number of iterations
%  [tau] : convergence tolerence
%
% OUTPUT:
%   x : nxm matrix of approximate solutions to Ax=b
% 
% Written by:
% -- 
% John L. Weatherwax                2006-12-19
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

if( nargin<3 | isempty(x) ) x=b; end
if( nargin<4 | isempty(M) ) M=1000; end
if( nargin<5 | isempty(tau) ) tau=sqrt(eps); end

% obtain the factors S and T in the splitting A = S - T: 
S=diag(diag(A)); T = S - A; 

divtol=100*max(abs(b-A*x));
for i=1:M
  x_new = S \ (T*x + b);           % the new value
  r=b-A*x_new;                     % the new residual
  dx=x_new - x;

  % convergence check
  if( all(abs(r)<tau) & all(abs(dx)<tau) ), return; end
  % divergence check
  if( max(abs(r))>divtol )
    disp('WARNING: Iterations diverging');
    return
  end
  
  x = x_new; % reassign the new x
end

if i==M, disp('WARNING: Maximum iterations exceeded'); end