function [x,n,nr,reshist]=it_gauss_seidel(b,A,N,tol)
% rescaling
D=diag(diag(A));
A1=inv(D)*A;
m=size(A,1);
b1=D\b;
nb=norm(b1);
% use 0 as initial guess
x=zeros(size(b));
for k=1:N
    for i=1:m
        s=0;
        for j=1:m
            if j~=i
                s=s+A1(i,j)*x(j);
            end
        end
        x(i) = b1(i)-s;
    end
    r=b1-A1*x;
    nr=norm(r);
    reshist(k)=nr/nb;
    if nr<tol*nb
        break;
    end
    x=x+r;
end
nr=norm(r)/nb;
n=k;
