  function [x,res]=sorRelax(A,x,b,maxit,omega,tol) 
% function [x,res]=sorRelax(A,x,b,maxit,omega,tol) 
% res stores consecutive residuals norms
 nn = size(A,1);
%%---------- compute residual
 r = b - A*x;
 res(1) = norm(r,1);
 tol = tol*res(1); 
%%---------- main loop
 for it=1:maxit
%%---------- coordinate loop
     for j=1:nn
%%---------- change to xj
       del = omega*r(j)/A(j,j); 
       x(j) = x(j) + del;
%%---------- update r accordingly
       r = r - del*A(:,j);
     end
     res(it) = norm(r,1) ;     
     if (res(it) < tol), break, end 
%% fprintf(1,'it = %4d  res-norm= %e \n',it,res(it))
 end
