function [x, iter, nf] = Lp_ls(A,b,sigma,p,lambda,tol,M,maxiter,x0,verbose)

% This code solves the following problem using a nonmonotone proximal
% gradient method
%
% min lambda*norm(A*x - b)^2 + sum(abs(x).^p)
%
% Input
% A, b        - function handle or m by n matrix and m vector from problem
% p           - between 0 and 1
% lambda      - penalty parameter
% tol         - tolerance for termination
% M           - integral parameter for nonmonotone line search
% maxiter     - maximum number of iterations
% x0          - a given initial point, not necessarily feasible
% verbose     - 1 if output is printed
%
% Output
% x           - approximate first-order stationary point
% iter        - total number of iterations (including inner iterations)
% nf          - total number of function evaluations
%
% Last updated:

if verbose == 1
  fprintf(' ***Inner loop: lambda = %g ***\n', lambda)
end

n = size(x0,1);
Lmin = 1;
Lmax = 1e8;
c = 1e-4;
tau = 2;
x = x0;

% Create function handle
if ~isa(A,'function_handle')
  A = @(x,mode) explicitMatrix(A,x,mode);
end

Axb = A(x,1)-b;
fx = norm(Axb)^2;
gx = A(Axb,2)*(2*lambda);
Fx = lambda*fx + sum((abs(x).^p));

Fl = -ones(M,1)*inf; 
Fl(1) = Fx; Fmax = Fx; 
L = 1; iter = 1; nf = 0;
err_lambda = inf;
diff_fx = inf;

logH = '%3s  %13s  %13s  %13s  %13s\n';
logB = '%3i  %13.6e  %13.6e  %13.6e  %13.6e\n';

% Print initial log.
if verbose == 1
  fprintf(logH,'itn','sm fval','L','|dx|','residual');
end

% Main loop
while iter <= maxiter && (err_lambda > sqrt(tol) || diff_fx > min(tol^2,1e-4) )
% terminate when iter is too large or err_lambda is small or change in function value very small
  
    Fold = Fx;
    while 1 == 1
        xnew = Lp_proj(x-gx/L,ones(n,1)/L,p);
        Axb = A(xnew,1)-b;
        fx = norm(Axb)^2;
        Fx = lambda*fx + sum((abs(xnew).^p));        
        nf = nf + 1;
        if (Fmax - Fx < c*norm(xnew-x)^2/2)
            L = tau*L;
        else
            break;
        end
    end
    gxnew = A(Axb,2)*(2*lambda);
    dx = xnew - x; 
    dg = gxnew - gx;
    
    % BB stepsize
    L = max(Lmin,min(Lmax,dx'*dg/norm(dx)^2));
	  Fl(mod(iter,M)+1) = Fx;
    Fmax = max(Fl);  
    x = xnew;
    gx = gxnew;
    
    % check for termination
    err_lambda = norm(x.*gx + p*(abs(x).^p), 'inf');
    diff_fx = abs(Fx - Fold)/max(1,abs(Fx));
    if mod(iter,10) == 0 && verbose == 1
      fprintf(logB, iter, Fx, L, norm(dx), max(norm(Axb)^2 - sigma^2, 0));
    end
    iter = iter + 1;
end

if verbose == 1
  fprintf(' *** Terminate inner loop: Itn = %3i  opt_cond = %13.6e\n', iter, err_lambda)
end


