function [x,fval,residual,iter,nf,lambda,time] = Lp2(A,b,sigma,p,tol,lambda0,eps0,rho,theta,M,maxiter,x0,verbose)

% This code solves the following problem using a penalty approach
%
% min sum(abs(x).^p) + lambda*norm(A*x - b)^2, 
%
% until max(norm(Axb)^2 - sigma^2, 0) = 0. The subproblem is solved using a nonmonotone proximal gradient
% scheme, with warmstart and penalty parameter updates
%
% Input
% A, b, sigma - m by n matrix, m vector and POSITIVE number from problem
% p           - between 0 and 1
% tol         - termination tolerance, for feasibility and first order
%               stationary condition
% lambda0     - initial penalty parameter
% eps0        - initial tolerance for solving the penalty subproblem
% rho         - rate of increase for lambda
% theta       - rate of decrease for eps
% M           - integral parameter for nonmonotone line search
% maxiter     - maximum number of iterations
% feas_pt     - a given feasible point
% x0          - a given initial point, not necessarily feasible
% verbose     - 1 if output is printed
%
% Output
% x           - approximate first-order stationary point
% fval        - sum(abs(x).^p)
% residual    - max(norm(A*x-b)^2 - sigma^2, 0)
% iter        - total number of iterations (including inner iterations)
% nf          - total number of function evaluations
% time        - time elapsed
%
% Last updated:

tstart = tic;
lambda = lambda0;
eps = eps0;
nf = 0; iter = 0;
x = x0;

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


while 1 == 1
    [x, iter_lambda, nf_lambda] = Lp_ls(A,b,sigma,p,lambda,eps,M,maxiter,x,verbose);
    iter = iter + iter_lambda; 
    nf = nf + nf_lambda;
    if max(norm(A(x,1)-b)^2 - sigma^2, 0) > tol || eps > tol*100 
        lambda = lambda*rho;
        eps = max(eps*theta,tol);
    else
        break;
    end
end

fval = sum(abs(x).^p);
residual = max(norm(A(x,1)-b)^2 - sigma^2, 0);
time = toc(tstart);
fprintf('\n\n --------------- The computational result ----------------\n\n')
fprintf(' Iter = %3.0d   nf=%d   fval = %7.3f   residual = %7.3f  nnz=%d   time = %.2f\n', ....
       iter, nf, fval, residual, nnz(x),time)
fprintf('\n ----------------------------------------------------------\n\n')

end

