function [u, ucache] = mmg_solver(A, b, lower, upper, u0, PP, maxnu, varargin)

m = length(PP)+1;
AA = cell(1,m);
AA{m} = A;
for k=m:-1:2
    AA{k-1} = PP{k-1}'*AA{k}*PP{k-1};
end

truncateCycle = 0;
if (length(varargin)>0)
    truncateCycle = varargin{1};
end

use_cache = (length(varargin)>1);
if (use_cache)
    ucache = zeros(size(A,1), maxnu+1);
    ucache(:,1) = u0;
end

u = u0;
for nu=1:maxnu
    r = b - A*u;
    ll = lower - u;
    uu = upper - u;
    if (truncateCycle==0)
        u = u + tmmg_step(r, ll, uu, AA, PP, 0);
    else
        u = u + tmmg_step(r, ll, uu, AA, PP, mod(nu-1,truncateCycle)==0);
    end
    if (use_cache)
        ucache(:,nu+1) = u;
    end
end
