function [X,iter,err] = GAUSSSEIDEL(A,B,epsilon,Niter,tol); 
% [X,iter,err] = GAUSSSEIDEL(A,B,epsilon,Niter,tol); 
[N,N] = size(A); TA = tril(A); FA = TA - A; 
ITA   = TA^(-1); AA = ITA*FA; 
X0    = 2*tol*(rand(N,1)-1/2); X = X0;  
err   = 2*epsilon; iter = 0;   
while (err > epsilon) && (iter < Niter)
   Xold = X; X = AA*X + ITA*B; 
   err  = norm(X-Xold); iter = iter + 1;
end
end