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.