function y = pgs_inc(A, r, lower, upper)

n = size(A, 1);
y = zeros(n,1);

for i = 1:n
    if (A(i,i)==0)
        continue;
    end
%    for j = find(A(i,:))
%        if (j<i)
%            r(i) = r(i) - A(i,j)*y(j);
%        end
%    end
    r(i) = r(i) - A(i,:)*y;
    y(i) = max([ min([ r(i)/A(i,i) upper(i) ]) lower(i)]);

%    y(i) = max([ min([ (r(i) - A(i,:)*y)/A(i,i) upper(i) ]) lower(i)]);
end
