function [x] = gauss_seidel(N, A,x,b);
%
%  one step of fixed point iteration for solving  A x = b
%
%  x^{j+1} = x^j  + M^{-1} (b - A x^j)
%
%  Gauss-Seidel:  M = D -L  (lower triangular part of A = D -L -U)
%
%  (D -L) x^{j+1} = (D -L) x^j  +(b - A x^j)
%                 = (D -L -A) x^j  +b
%                 = U x^j +b
%
% ==>  D x^{j+1}  = L x^{j+1}  +U x^j  +b



%%%

%  %  direct implementation
%  %  only for N small, e.g. N < 10
%
%  for i = 1 : N
%  	x(i) = b(i);
%  	for k = 1 : i-1
%  		x(i) = x(i) -A(i,k)*x(k);	% L x^{j+1}
%  		end
%  	for k = i+1 : N
%  		x(i) = x(i) -A(i,k)*x(k);	% U x^j
%  		end
%  	x(i) = x(i) / A(i,i);
%  	end

%%%

%  %  replace slow for loops with colon operator :
%
%  for i = 1 : N
%  	x(i) = ( b(i) -A(i,1:i-1)*x(1:i-1) -A(i,i+1:N)*x(i+1:N) ) / A(i,i);
%  	end

%%%

%  use only matrix computations

L = -tril(A,1);	% strict lower triangular part of A
U = -triu(A,1);	% strict upper triangular part of A
D = A +L +U;	% diagonal of A

x = (D-L) \ (U*x +b);
