function x=symmetricGSeidel(A,f,x,M,tau);
% GSEIDEL: Solves Ax=f using the symmetric Gauss-Seidel iteration.
% USAGE:
%   x=symmetricGSeidel(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);

% obtain the coefficient matrix (and right hand side) when sweeping in 
%   ASSENDING order of the unknowns: 
PG_A = (D-L) \ U; 
f_A  = (D-L) \ f; 

% obtain the coefficient matrix (and right hand side) when sweeping in 
%   DECENDING order of the unknowns: 
PG_D = (D-U) \ L; 
f_D  = (D-U) \ f; 

divtol=100*max(abs(f-A*x));
for i=1:M
  if( mod(i,2)~=0 ) % sweep in ascending order (D-L) x^{n+1} = U x^n + f
    x_new = PG_A * x + f_A; 
  else              % sweep in descending order (D-U) x^{n+1} = L x^n + f 
    x_new = PG_D * x + f_D; 
  end
  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