function [y,iter,L,stopcrit]=RFISTA(F,PAR)
%% unconstrained minimization

%%% solver specific parameters
PAR.A0=0;
PAR.tau0=1;

if isfield(PAR,'maxiter')
    maxj = PAR.maxiter;
else
    maxj = 1e6;
end

tildemu=PAR.mu;
tau=PAR.tau0;
L=PAR.Munderbar;

%%%% data
x=PAR.x0;
y=PAR.x0;
A=PAR.A0;

j = 0;       % major R-FISTA iterations
iter=0;      % prox evaluations
stopcrit = inf;
FuncValue = inf;

%%%%%%%%the big main outer loop starts here
while stopcrit > PAR.epsilon  && j <= maxj

    j = j + 1;
    yold = y;

    %%% step 1 with Lj == Lj-1
    a=(tau+sqrt(tau^2+4*tau*A*L))/(2*L);
    xtilde=((A*yold)+(a*x))/(A+a);
    storegradient=F.grad(xtilde);
    y=F.prox(xtilde-(1/L)*storegradient,1/L);
    iter=iter+1;

    FuncValuextilde=F.fs(xtilde);
    FuncValue=F.fs(y);
    NormDiff=PAR.norm_fn(y-xtilde)^2;
    LinearTerm=PAR.prod_fn(y-xtilde,storegradient);

    Ftest =  FuncValuextilde + LinearTerm + ...
        ((1-PAR.chi)*L*NormDiff)/(4)+10^-6 - FuncValue;

    while Ftest < 0
        L=PAR.beta*L;

        a=(tau+sqrt(tau^2+4*tau*A*L))/(2*L);
        xtilde=((A*yold)+(a*x))/(A+a);
        storegradient=F.grad(xtilde);
        y=F.prox(xtilde-(1/L)*storegradient,1/L);
        iter=iter+1;

        FuncValuextilde=F.fs(xtilde);
        FuncValue=F.fs(y);
        NormDiff=PAR.norm_fn(y-xtilde)^2;
        LinearTerm=PAR.prod_fn(y-xtilde,storegradient);

        Ftest =  FuncValuextilde + LinearTerm + ...
            ((1-PAR.chi)*L*NormDiff)/(4)+10^-6 - FuncValue;
    end

    if j==1
        PAR.mu=PAR.muEst;
        %PAR.mu=(4*(FuncValue-FuncValuextilde-LinearTerm))/((1-PAR.chi)*NormDiff);
        tildemu=PAR.mu;
        % keyboard
    end

    %%% step 2
    s=L*(xtilde-y);
    v=F.grad(y)-storegradient+s;
    A=A+a;
    oldtau=tau;
    tau=oldtau+(0.5*tildemu*a);
    x=(1/tau)*((0.5*tildemu*a*y)+(oldtau*x)-(a*s));

    %%% step 3 restart
    if (PAR.norm_fn(y-PAR.x0))^2<PAR.chi*(A*L*NormDiff)
        PAR.Mbar=L;
        PAR.Munderbar=max(0.4*PAR.Mbar,PAR.Mbar0);
        L=PAR.Munderbar;
        PAR.mu=0.5*tildemu;
        tildemu=PAR.mu;

        PAR.x0=y;
        x=PAR.x0;
        y=PAR.x0;

        PAR.A0=0;
        A=PAR.A0;

        PAR.tau0=1;
        tau=PAR.tau0;
    end

    stopcrit=PAR.norm_fn(v);
end
% keyboard
% PAR.Munderbar=L/2;
% PAR.Munderbar=10;
PAR.Mbar0=PAR.LipEst/4;
PAR.Munderbar=PAR.Mbar0;

% keyboard
end
