function [fval,z,x,y,iter] = DR_proj_raw(A,b,k,tol,maxiter,freq)

% This computes the intersection of 
%
% ||u||_0 <= k, ||u||_infty <= 1e6 and Au = b, A has full row rank, m < n
%
% Using DR splitting to minimize the sum of the indicator function of C and
% the sparsity constraint, where C is the affine space.

n = size(A,2);

% initialization
AAT = A*A';
x = zeros(n,1);
y = x;
z = x;

iter = 0;
while iter <= maxiter
  
  xold = x;
  yold = y;
  zold = z;
  
  % update y: proj onto Ax = b
  y = x + A'*(AAT\(b - A*x)); % projecting x onto A*x = b

  % 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));
  
  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 mod(iter,freq) == 0
    projz = z + A'*(AAT\(b - A*z)); % Just for checking fval, not necessary for convergence
    fval = 0.5*norm(z - projz)^2;
    fprintf(' iter %d fval %g err %g normx %g\n', iter, fval, err, norm(x))
  end
  
end
projz = z + A'*(AAT\(b - A*z)); % Just for checking fval, not necessary for convergence
fval = 0.5*norm(z - projz)^2;
fprintf(' Termination: ')
fprintf(' iter %d fval %g err %g\n', iter, fval, err)