function s=DSS(V,U,r);
%V is k x n and U is l x n have orthonormal rows.
%VV^T = I_k and UU^T=I_l
%k<r, we will choose k columns.
%s are the scaling weights for each column (the non-zeros).
%The sampled, rescaled columns of U and V satisfy Lemma 3.5 in the paper: 
%http://arxiv.org/abs/1103.0995
%If U=V, you obtain the single set sparsification result of Batson, Spielman, 
%srivastava, 2009.
%
%Copyright Malik Magdon-Ismail


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


s=zeros(n,1);
Atau=zeros(k,k);
Btau=zeros(l,l);
delL=1;
delU=(1+sqrt(l/r))/(1-sqrt(k/r));
Ik=eye(k);
Il=eye(l);

for tau=0:r-1;
	Ltau=tau-sqrt(r*k);
	Utau=delU*(tau+sqrt(r*l));
	ALinv=inv(Atau-(Ltau+delL)*Ik);
	ALinv2=ALinv*ALinv;
	lamAtau=svd(Atau);
	phiL1=sum(1./(lamAtau-Ltau-delL));
	phiL2=sum(1./(lamAtau-Ltau));
	BUinv=inv((Utau+delU)*Il-Btau);
	BUinv2=BUinv*BUinv;
	lamBtau=svd(Btau);
	phiU1=sum(1./(Utau-lamBtau));
	phiU2=sum(1./(Utau+delU-lamBtau));
	for i=1:n
		vi=V(:,i);
		ui=U(:,i);
		Lower=vi'*ALinv2*vi/(phiL1-phiL2)-vi'*ALinv*vi;
		Upper=ui'*BUinv2*ui/(phiU1-phiU2)+ui'*BUinv*ui;
		if(Upper<=Lower)
			t=2/(Upper+Lower);
			s(i)=s(i)+t;
			Atau=Atau+t*vi*vi';
			Btau=Btau+t*ui*ui';
			break;
		end
	end
end

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

