function [xt,dent,ft,void,cft] = ...
    armijo_QNP(x,step,gn,den,g0,f0,eta,TF,bg,epsilonk)
%ARMIJO_QNP   Armijo rule along the projection arc.
%   ARMIJO_QNP computes the new QNP iterate using the Armijo rule along
%   the projection arc for the step-length computation.
%
%   [XT]=ARMIJO_QNP(X,STEP,GN,DEN,G0,F0,ETA,TF,BG,EPSILONK) returns 
%   the new QNP iterate XT. X is the current iterate. STEP is the search
%   direction. GN is the data image. DEN is A*X+BG. F0 is the cost function
%   at the current iterate and G0 is its gradient. ETA is the 
%   regularization parameter. TF is the fft2 of the psf image. BG is the 
%   background value. EPSILONK is the parameter for the set Ik. 
%
%   [XT,DENT]=ARMIJO_QNP(X,STEP,GN,DEN,G0,F0,ETA,TF,BG,EPSILONK) 
%   also returns the value A*XT+BG at the new iterate.
%
%   [XT,DENT,FT]=ARMIJO_QNP(X,STEP,GN,DEN,G0,F0,ETA,TF,BG,EPSILONK) 
%   also returns the objective function value FT at the new iterate.
%
%   [XT,DENT,FT,VOID]= ...
%   ARMIJO_QNP(X,STEP,GN,DEN,G0,F0,ETA,TF,BG,EPSILONK)
%   also returns a parameter VOID:
%    0 the line search terminates successfully
%    1 the line search terminates after 50 backtracks without determining 
%    a new iterate.
%
%   [XT,RESIDT,FT,VOID,CFT]=ARMIJO_QNP(X,STEP,GN,G0,F0,ETA,EPSILONK)
%   also returns the number CFT of performed FFTs.

%   D. P. Bertsekas. On the Goldstein-Levitin-Poljak Gradient Projection 
%   Method. IEEE Transcations on Automatic Control, 21(2):174–184, 1976

cft = 0; % number of ffts
maxarm = 50; % maximum number of backtracks
iarm = 1;
void = 0;
lambda = 1; 
alpha = 1e-4;
beta = 0.5;

xt = max(x+lambda*step,0);      
dent = real(ifft2(TF.*fft2(xt)))+bg; cft = cft+2;                   
ft = sum(sum(gn.*log(gn./dent)+dent-gn))+eta*0.5*norm(xt,'fro')^2;
 
Ik = (x<=epsilonk & g0>0);
stept = lambda*step;
stept(Ik) = -(x(Ik)-xt(Ik));

while (ft >= f0 + alpha * g0(:)'* stept(:));
    lambda = lambda*beta;
    xt = max(x+lambda*step,0); 
    dent = real(ifft2(TF.*fft2(xt)))+bg; cft=cft+2;
    ft = sum(sum(gn.*log(gn./dent)+dent-gn))+eta*0.5*norm(xt,'fro')^2;
    stept = lambda*step;
    stept(Ik) = -(x(Ik)-xt(Ik));
    iarm = iarm+1;
    if iarm > maxarm
        void = 1;
        return;
    end
end  % end line search


