function s=DSF(V,A,r)
%V is k x n with orthonormal rows A is l x n .
%VV^T = I_k
%k<r, we will choose r columns.
%s are the scaling weights for each column (the non-zeros).
%The sampled, rescaled columns of V and A satisfy Lemma 3.6 in the paper: 
%http://arxiv.org/abs/1103.0995
%
%Copyright Malik Magdon-Ismail

[k,n]=size(V);
[l,m]=size(A);
if(m~=n)
	s='dimensions do not match'
	return
end

normA=sum(A.^2);

jpos=find(normA>0);

s=zeros(n,1);
Atau=zeros(k,k);
delL=1;
delU=sum(normA)/(1-sqrt(k/r));

Ik=eye(k);

for tau=0:r-1;
	Ltau=tau-sqrt(r*k);
	Utau=delU*tau;
	ALinv=pinv(Atau-(Ltau+delL)*Ik);
	ALinv2=ALinv*ALinv;
	lamAtau=svd(Atau);
	phiL1=sum(1./(lamAtau-Ltau-delL));
	phiL2=sum(1./(lamAtau-Ltau));
    	LA=zeros(1,n); UA=zeros(1,n);
	for i=1:n
		vi=V(:,i);
		Lower=vi'*ALinv2*vi/(phiL1-phiL2)-vi'*ALinv*vi;
        	LA(i)=Lower;
        	Upper=normA(i)/delU;
        	UA(i)=Upper;
    	end
	candidates0 = find(UA<=LA);
	candidates = intersect(jpos,candidates0);	
	nnz_s = find(s'>0);
	Vsel = V(:,nnz_s);
	if(~isempty(candidates))
		index_diff = setdiff(candidates,nnz_s);
        	if (isempty(index_diff))
			Vsel = V(:,candidates);
            		colnorm = sqrt(sum(Vsel.^2,1));
            		[val,max_idc] = max(colnorm);
            		col_sel = candidates(max_idc);
        	else
			Vsel = V(:,index_diff);
            		colnorm = sqrt(sum(Vsel.^2,1));
            		[val,max_idc] = max(colnorm);
            		col_sel = index_diff(max_idc);
        	end

        	t = 2/(UA(col_sel)+LA(col_sel));
 		s(col_sel) = s(col_sel) + t;
        	vj = V(:,col_sel);
        	Atau = Atau + (t * vj * vj');
    	else
		disp('Something is wrong: empty candidate set'); 
    	end
%	[tau,r,max(s),sum(s.*normA'),Utau+delU,svds(Atau,1,0),Ltau+1]
%	pause
    
%         if(Upper<=Lower)
%             
% 			t=2/(Upper+Lower);
% 			s(i)=s(i)+t;
% 			Atau=Atau+t*vi*vi';
% 			break;
%         end 
end

s=s*(1-sqrt(k/r))/r;
