Methods for simultaneous diagonalization?

I have two matrices that commute (within machine precision), both of which have eigenvalues that may be either degenerate or very closely spaced. Is there an existing implementation of a simultaneous diagonalization routine in Julia that gives a stable way of computing the eigenvectors? I’m afraid that the (near-)degeneracy will cause serious issues…

To be a bit more specific, I am trying to compute the Bloch-Messiah decomposition of a unitary matrix, which involves putting an antisymmetric matrix into canonical form (block diagonal where each 2x2 block is of the form [0.0, value; -value, 0.0]). I am trying to see if simultaneous diagonalization can be used here to stably generate the required decomposition.

(See also GitHub - JuliaLinearAlgebra/SkewLinearAlgebra.jl: Julia Package for skew-symmetric matrices for optimized eigensolvers of anti-symmetric matrices. Given the eigen-decomposition you could easily compute the block-diagonal form.)