function x=wJacobi(A,f,omega,x,M,tau);
% WJACOBI: Solves Ax=f using weighted Jacobi iterations
% USAGE:
%   x=WJACOBI(A,f,x,M,tau,lambda);
%
% INPUTS:
%   A : nxn matrix
%   f : nxm matrix
%  [omega]: relaxation parameter. Default: (2/3) 
%           (omega=1 gives the original Jacobi method)
%  [x] : nxm matrix of starting values. Default: f
%  [M] : maximum number of iterations. Default: 1000
%  [tau] : convergence tolerence. Default: sqrt(eps)
%
% OUTPUT:
%   x : nxm matrix of approximate solutions to Ax=f
% 
% 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(omega) ) omega=2/3; end
if( nargin<4 | isempty(x) ) x=f; end
if( nargin<5 | isempty(M) ) M=1000; end
if( nargin<6 | isempty(tau) ) tau=sqrt(eps); end

% obtain the factors D, L and U in the splitting A = D - L - U:
D=diag(diag(A)); L = -tril(A,-1); U = -triu(A,1);

% compute the Jacobi iteration matrix: 
PJ = D \ (L+U); 

% compute the modified right hand side:
omegaF = omega * (D \ f); 

% compute the weighted Jacobi iteration matrix: 
wPJ = (1-omega)*eye(size(A)) + omega*PJ; 

divtol=100*max(abs(f-A*x));
for i=1:M
  x_new = wPJ * x + omegaF;        % the new value
  r=f-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