function xnew = pgs(A, b, lower, upper, xold)

r = b - A*xold;
xnew = xold + pgs_inc(A, r, lower-xold, upper-xold);

