  function [x,nrmHist] = sor (A, b,x, nsteps, om, tol)
 D = diag(diag(A));
 E = -tril(A,-1);
 F = -triu(A,1);
 L = D-om*E;
 U = om*F + (1-om)*D;
%%-------------------- record res. norm
 ro = norm(b - A*x);
 nrmHist(1) = ro;
 tol_rel  = tol* ro; 
%%-------------------- iterate 
for i=1: nsteps
  x = L \ (U*x + om*b);
  ro = norm(b - A*x); 
  nrmHist(i+1) = ro;
  if (ro < tol_rel)
    break;
  end
end

