function [x,nit,err]=SOR(A,b,theta)
x0=b;
eps0=1.e-8;
epsb=eps0*norm(b);
err=1;
n=length(b);
x=zeros(n,1);
for nit=1:1000
    for i=1:n
        j=1:i-1; k=i+1:n;
        x(i)=theta*((b(i)-A(i,j)*x(j)-A(i,k)*x0(k))/A(i,i))+(1-theta)*x0(i);
    end
    r=b-A*x;
    inc=x-x0;
    if norm(r)<epsb && norm(inc)<epsb
        err=0;
        return
    end
    x0=x;
end