function [L, D, R] = odo(U)

%odo=orthogonal diagonal orthogonal
%Given a  unitary matrix U,
%it finds special (det=1) orthogonal matrices L and R
%and a diagonal unitary matrix D
%such that L*D*R'=U

dim_U= rows(U);
if(norm(U*U'-eye(dim_U))>1e-8)
	error("input matrix for odo decomposition is not unitary");
end


X = (U + conj(U))/2;%IMP: this is NOT equal to (U+U')/2
Y = (U - conj(U))/(2*i);
% U = X + i*Y; X, Y both real

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

D = DX + i*DY;