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.