function [ x,iter,talg ] = IRL1e3( A,b,epsilon,lambda, L, tol,maxiter,re_freq)
% This function computes a stationary point for the following model
% 
% min 0.5*||Ax - b||^2 + lambda*sum((abs(x)/epsilon + 1))
% 
% It uses IRL1e3
%
% Input
%
% A       - m by n matrix (m << n)
% b       - m by 1 vector measurement
% lambda  - real number > 0
% epsilon - real number > 0
% L       - user supply value of max(eig(A*A'))
% tol     - tolerance
% maxiter - max number of iteration
% re_freq - iteration beyond which theta is not updated
% 
% Output
%
% x       - approximate stationary point
% iter    - number of iterations
% talg    - cputime

tstart_alg = tic;
[~,n] = size(A);
x = zeros(n,1);
z = x;
iter = 1;
ell = lambda/epsilon^2;
%generate theta_0
theta = 1;
for count = 1:6
    theta = 2/(1 + sqrt(1+4/theta^2));     
end


%updates in IRL1e3
while iter < maxiter
    %update s
    s = lambda./(abs(x) + epsilon);
    
    %update y
    y = (1-theta)*x + theta*z;
    
    % update z
    G = A'*(A*y-b);
    tmp = z - G/theta/L;
    z = sign(tmp).*max(abs(tmp) - s/theta/L,0);
    
    %update x
    xold = x;
    tmp = y - G/L;
    x = sign(tmp).*max(abs(tmp) - s/L,0);
    
    %check for termination
    if (2*L*norm(x - y) + ell*norm(x - xold))/max(1,norm(x)) < tol
        break
    end
    
    %update theta
    if iter <= re_freq
        theta = 2/(1 + sqrt(1+4/theta^2));        
    end
    
    iter = iter + 1;    
end

talg = toc(tstart_alg);