function [x, iter]= sor ( A, b, x0, nmax, toll, omega)
[~,n]=size(A);
iter = 0; r = b - A * x0; r0 = norm(r); err = r0; x=x0; 
D=diag(diag(A)); S=tril(A,-1)+(1/omega)*D; T=-triu(A,1)+(1-omega)/omega*D;  
while err > toll && iter < nmax
iter = iter + 1;
x = S\(T*x+b);
err=norm(b-A*x)/r0;
end