function [U,V,D] = Chenmh_CSD(Q,r);
% %% Q should be column-orthonormal
[m,n] = size(Q); % m>=n
D = zeros(m,n);
[U1,C,V] = svd(Q(1:r,1:n));
D(1:r,1:n) = C;

S = diag(diag(eye(n) - C'*C).^0.5);
D(m-n+1:m,1:n) = S;
D(1:r,1:n) = C;
Temp1 = Q(r+1:m,1:n)*V*D(r+1:m,1:n)';
[U2,R1] = qr(Temp1(:,end:-1:1)); 
ms1 = (diag(R1))<-eps; ms1 = diag(1-2.*ms1); 
U2 = U2*ms1; U2 = U2(:,end:-1:1);

U = [U1 zeros(r,m-r);zeros(m-r,r) U2];