function [fval,z,x,y,iter] = DR_ls(A,b,k,tol,maxiter,ratio)
% This computes the constrained least squares problem
% min 0.5*||Au - b||^2  subject to ||u||_0 <= k, ||u||_infty <= 1e6

% using DR splitting algorithm 

[m,n] = size(A);


% initialization
AAT = A*A';
nmA = norm(AAT);
Ab = A'*b;
x = zeros(n,1);
y = x;
z = x;
gamma0 = (sqrt(1.5) - 1)/nmA; % has convergence guarantee when gamma < gamma0
gamma = gamma0*ratio; % heuristic

iter = 0;
while iter <= maxiter

  xold = x;
  yold = y;
  zold = z;
  
  % update y: 
  tmp = x + gamma*Ab;
  y = tmp - gamma*(A'*((eye(m) + gamma*AAT)\(A*tmp)));
%   fprintf(' Check error %g\n', norm(y + gamma*(A'*(A*y)) - tmp))
  
  % update z
  tmp = 2*y - x;
  tmp0 = max(min(tmp,1e6),-1e6);
  tmp2 = tmp.^2 - (tmp0 - tmp).^2;
  [~,I] = sort(tmp2,'descend');
  z = zeros(n,1);
  z(I(1:k)) = tmp0(I(1:k));
  
  % update x
  x = x + z - y;
  
  iter = iter + 1;
  err = max([norm(x - xold), norm(y - yold), norm(z - zold)]);
  if err/max([norm(xold),norm(yold),norm(zold),1]) < tol
    break
  end

  if (norm(y - yold) >= 1e3/iter || norm(y) > 1e10) && gamma > gamma0 % decrease gamma to guarantee convergence
    gamma = max(gamma/2,gamma0*0.9999);
  end
  
%   if mod(iter,freq) == 0
%       fval = 0.5*norm(A*z - b)^2;
%       fprintf('iter %d fval %g  err %g  gamma %g normx %g erry %g\n', iter, fval, err, gamma, norm(x),norm(y - yold))
%   end
      
end

fval = 0.5*norm(A*z - b)^2;
fprintf(' Termination: ')
fprintf(' iter %d fval %g err %g\n', iter, fval, err)
  
  
  
end

