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

use_cache = (length(varargin)>0);
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;
    u = u + mmg_step(r, ll, uu, AA, PP);
    if (use_cache)
        ucache(:,nu+1) = u;
    end
end
