function x=redBlackGSeidel(A,f,x,M,tau);
% GSEIDEL: Solves Ax=f using Red-Black Gauss-Seidel iteration.
% USAGE:
%   x=redBlackGSeidel(A,f,x,M,tau);
%
% INPUTS:
%   A : nxn matrix
%   f : nxm matrix
%  [x] : nxm matrix of starting values. Default: f
%  [M] : maximum number of iterations
%  [tau] : convergence tolerence
%
% 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(x), x=f; end
if nargin<4 | isempty(M), M=1000; end
if nargin<5 | 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); I = eye(size(A)); Z = zeros(size(f)); 

% obtain the coefficient matrix when processing the
%   EVEN unknowns (the identity and zero avoid touching the ODD unknowns): 
PG_E = (D-L) \ U; PG_E(1:2:end,:) = I(1:2:end,:); 
fG_E = (D-L) \ f; fG_E(1:2:end,:) = Z(1:2:end,:); 

% obtain the coefficient matrix when processing the
%   ODD unknowns (the identity and zero avoid touching the EVEN unknowns): 
PG_O = (D-L) \ U; PG_O(2:2:end,:) = I(2:2:end,:); 
fG_O = (D-L) \ f; fG_O(2:2:end,:) = Z(2:2:end,:); 

divtol=100*max(abs(f-A*x));
for i=1:M
  % first update the even components
  x_new = PG_E * x + fG_E;
  % second update the odd components
  x_new = PG_O * x_new + fG_O; 

  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; 
end

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