function x=gseidel(A,b,x,M,tau,lambda);
% GSEIDEL: Solves Ax=b using Gauss-Seidel iteration.
% USAGE:
%   x=gseidel(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 = tril(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;                     % 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