I would like to diagonalize a very large matrix that has the same property as quaternion matrices (in a sense that the matrix can be written as a linear combination of quaternions): $$ H = \left(\begin{array}{cc} H_{11} & H_{12} \\ -H^\ast_{12} & H^\ast_{11} \end{array}\right), \quad U^\dagger HU = \epsilon $$ where $\epsilon$ is a diagonal matrix, $U$ is a unitary matrix and $U$ has the same symmetry $$ U = \left(\begin{array}{cc} U_{11} & U_{12} \\ -U^\ast_{12} & U^\ast_{11} \end{array}\right) $$ Note that there are $dim(\epsilon)/2$ sets of identical eivenvalues.
Of course it is trivial to diagonalize $H$ in a brute-force way to get right eivenvalues; however, for some reasons I need a symmetry-adapted set of eivenvectors as shown above. Is there a way to transform the matrix obtained by brute-force diagonalization to an symmetry-adapted matrix above? Although I am aware of some papers such as (J. Math. Phys. 46, 052106 (2005); http://dx.doi.org/10.1063/1.1896386), which suggest $$ \left(\begin{array}{cc} V_{11} & V_{12} \\ V_{21} & V_{11} \end{array}\right)\to \left(\begin{array}{cc} U_{11} & U_{12} \\ -U^\ast_{12} & U^\ast_{11} \end{array}\right) $$ where $V^\dagger H V = \epsilon$ and $$ U_{11} = \frac{1}{2}(V_{11}+V_{22}^\ast),\quad U_{12}=\frac{1}{2}(V_{12}-V_{21}^\ast) $$ it does not seem to preserve orthogonality of the eigenvectors.
