A symplectic matrix is a $2n\times2n$ matrix $S$ with real entries that satisfies the condition
$$ S^T \Omega S = \Omega $$ where $\Omega$ is the symplectic form, typically chosen to be $\Omega=\left(\begin{smallmatrix}0 & I_N \\ -I_N & 0\end{smallmatrix}\right)$. Sympletic matrices form the symplectic group $Sp(2n,\mathbb{R})$. Any symplectic matrix S can be decomposed as a product of three matrices as
\begin{equation} S = O\begin{pmatrix}D & 0 \\ 0 & D^{-1}\end{pmatrix}O' \quad \quad \forall S \in Sp(2n,\mathbb{R}), \end{equation} where $O, O'$ are orthogonal and symplectic - $\operatorname{Sp}(2n,\mathbb{R})\cap \operatorname{O}(2n)$; $D$ is positive definite and diagonal. The form of a matrix that is both symplectic and orthogonal can be given in block form as $O=\left(\begin{smallmatrix}X & Y \\ -Y & X\end{smallmatrix}\right)$, where $XX^T+YY^T=I_N$ and $XY^T-YX^T=0$. The decomposition above is known as Euler decomposition or alternatively as Bloch-Messiah decomposition.
How can I find the matrices in the decomposition for a given symplectic matrix?
Apparently, the decomposition is closely related to the singular value decomposition and I think the elements of the matrices $D$ and $D^{-1}$ coincide with the singular values of $S$. Also, I have the impression that the case where it can be assumed that $S$ is also symmetric is easier. Any help, tips or pointers would be much appreciated!