function [u,k] = Jacobi_2blocks(A,b,u0,I1,I2,prec)

D1 = inv(full(A(I1,I1)));
D2 = inv(full(A(I2,I2)));

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

while normr/normr0 > prec
    u(I1) = u(I1) + D1*r(I1);
    u(I2) = u(I2) + D2*r(I2);
    normr_prev = normr;
    r = b - A*u;
    normr = norm(r);
%     normr/normr_prev, normr/normr0, pause
    k = k+1;
end