function [L, DX, DY, R] = simul_real_svd(X,Y)

%This funtion does a
%simultaneous real singular value decomposition of 2 matrices.
%Given two real matrices X,Y of the same dimension,
%it finds special (det=1) orthogonal matrices L and R
%and real diagonal matrices DX, DY
%such that L*DX*R'=X and L*DY*R'=Y

%Reference: C.Eckart, G.Young, Bull. Am. Math. Soc. 45(1939)118


%LX'*X*RX = [DX00, Z01]
%           [Z10 , Z11]

%LX'*Y*RX = G = [G00, Z01]
%               [Z10, G11]
 
%[P00', Z10' ]*[DX00, Z01]*[P00, Z01 ] = DX
%[Z01', LG11'] [Z10 , Z11] [Z10, RG11]

%[P00', Z10' ]*[G00 , Z01]*[P00, Z01 ] = [DG00, Z01] = DY
%[Z01', LG11'] [Z10 , G11] [Z10, RG11]   [Z10, DG11]

dim_U= rows(X);
if(rows(X)!=rows(Y) | columns(X)!=columns(Y))
	error("X and Y don't have the same dimensions");
endif

[LX,DX,RX]=svd(X); %LX*DX*RX'=X

%next find rank of X
tol = dim_U*DX(1,1)*eps;
rank_X = 1;
j=2;
done=false;
while(j<dim_U+1 & !done)
	if(DX(j,j)>tol)
		rank_X = rank_X +1;
	else
		done=true;
	end
	j=j+1;
end

if(rank_X<dim_U)
	DX00 = DX(1:rank_X, 1:rank_X);
	Z01 = zeros(rank_X, dim_U-rank_X);
	Z10 = zeros(dim_U-rank_X, rank_X);
	%Z11 = zeros(dim_U-rank_X, dim_U-rank_X);

	G = LX'*Y*RX;
	G00 = G(1:rank_X, 1:rank_X);
	G11 = G(rank_X+1:dim_U, rank_X+1:dim_U);

	[ P00 ,  qDs ] = rjd([DX00, G00], 1e-10);
	DG00 = qDs(1:rank_X, rank_X+1: 2*rank_X);
	%P00*DX00*P00'=DX00
	%P00*DG00*P00'=G00

	[LG11,DG11,RG11]=svd(G11); %LG11*DG11*RG11' = G11
	DY = [DG00, Z01; Z10, DG11];

	R = RX*[P00,Z01;Z10,RG11];
	L = LX*[P00,Z01;Z10,LG11];
	%L*DX*R'=X
	%L*DY*R'=Y
else
	DX00 = DX;
	G00 = LX'*Y*RX;
	[ P00 ,  qDs ] = rjd([DX00,G00], 1e-12);
	DY = qDs(1:dim_U, dim_U+1: 2*dim_U);
	R = RX*P00;
	L = LX*P00;
end

if(det(L)<0)
	L(:,1) = -L(:,1);
	DX(1,1) = -DX(1,1);
	DY(1,1) = -DY(1,1);
end
if(det(R)<0)
	R(:,1) = -R(:,1);
	DX(1,1) = -DX(1,1);
	DY(1,1) = -DY(1,1);
end