function [u,k] = Jacobi_Nblocks(A,b,u0,I,prec)

N = length(I);
for i=1:N
    D{i} = inv(full(A(I{i},I{i})));
end

u = u0;
r = b-A*u;
normr = norm(r);
normr0 = normr;
k = 0;

while normr/normr0 > prec
    for i=1:N
        u(I{i}) = u(I{i}) + D{i}*r(I{i});
    end
    normr_prev = normr;
    r = b - A*u;
    normr = norm(r);
%     normr/normr_prev, normr/normr0, pause
    k = k+1;
end