function [x, xcache] = pgs_solver(A, b, lower, upper, x0, maxnu, varargin)

use_xcache = (length(varargin)>0);
if (use_xcache)
    xcache = zeros(size(A,1), maxnu+1);
    xcache(:,1) = x0;
end

xold = x0;
for nu=1:maxnu
    x = pgs(A, b, lower, upper, xold);
    xold = x;
    if (use_xcache)
        xcache(:,nu+1) = x;
    end
end
