function [time,NormResidualz,NormResidualR]=APFIAL(F,PAR,z,p)

if ~isfield(PAR,'verbose')
    PAR.verbose = true;
end
Residualz=F.Residual(z);
NormResidualz=norm(Residualz);
NormResidualR=inf;
totalIter=0;
k=0;

if PAR.verbose
    fprintf('\nAPFIAL\n')
    fprintf('%6s %10s %12s %14s %14s %14s %14s\n', ...
        'outer','inner','total','primal_res','inner_res','epsilon','c')
end

tic
while NormResidualz>PAR.epsilonhat || NormResidualR>PAR.rhohat
    %%% Index outer iteration
    k=k+1;
    PAR.k = k;

    %%% Tolerance for RFISTA method
    if NormResidualz>PAR.epsilonhat
        PAR.epsilon=max(PAR.tildepsilon,PAR.epsilonhat);
    else
        PAR.epsilon=PAR.epsilonhat;
    end

    %%% Update AL oracle for this penalty/multiplier pair.
    [fs,grad] = makeAugLagOracle(F,PAR,p,PAR.c);
    F.fs = fs;
    F.grad = grad;

    %%% Update initial point fed into RFISTA
    PAR.x0=z;
    %keyboard
    [z,iter,PAR.Munderbar,NormResidualR]=RFISTA(F,PAR);
    %keyboard
    totalIter=iter+totalIter;

    %%% Update residuals
    Residualz=F.Residual(z);
    NormResidualz=norm(Residualz);

    %%% Lagrange multiplier update
    p=p+(PAR.c*Residualz);

    if PAR.verbose
        fprintf('%6i %10i %12i %14.6e %14.6e %14.6e %14.6e\n', ...
            k, iter, totalIter, NormResidualz, NormResidualR, PAR.epsilon, PAR.c)
    end

    %%% Update tolerance and penalty parameter
    PAR.c=PAR.alpha*PAR.c;
    PAR.tildepsilon=PAR.tildepsilon/PAR.omega;
end
time=toc;
outerIter=k;

if PAR.verbose
    fprintf('Finished APFIAL: outer=%d, total inner=%d, time=%.2f seconds.\n', ...
        outerIter, totalIter, time)
end

end

function [fs,grad] = makeAugLagOracle(F,PAR,p,c)
finit = F.finit;
gradinit = F.gradinit;
Residual = F.Residual;
Adjoint = F.Adjoint;
prod_fn = PAR.prod_fn;
norm_fn = PAR.norm_fn;

last_x = [];
last_r = [];

fs = @valueOracle;
grad = @gradOracle;

    function value = valueOracle(x)
        r = residualAt(x);
        value = finit(x) + prod_fn(p,r) + 0.5*c*norm_fn(r)^2;
    end

    function g = gradOracle(x)
        r = residualAt(x);
        g = gradinit(x) + Adjoint(p + c*r);
    end

    function r = residualAt(x)
        if isempty(last_x) || ~isequal(x,last_x)
            last_x = x;
            last_r = Residual(x);
        end
        r = last_r;
    end
end
