function M = sparsematrix(M,s)
%% projection onto D = {B : max_{j=1,..,n} ||B_j||_0 <= s}, where B_j is the j_th column of B, B is an n by n matrix.

% Modified on Feb 19: column sparse
[n,~] = size(M);
[~,ind] = sort(abs(M),'ascend');
tmp = ind(1:n-s,:) + repmat(0:n:n*(n-1),n-s,1);
M(tmp(:))=0;
