function [x,iter]=jor(A,b,x0,nmax,toll,omega)
iter=0; r=b-A*x0; r0=norm(r); err=r0; x=x0; d=diag(A); A=A-diag(d);
while err>toll && iter<nmax
iter = iter + 1;
x=(1-omega)*x-omega*(A*x-b)./d;
r=b-(A+diag(d))*x; err=norm(r)/r0;
end