  It is well known that a set of non-defect matrices can be simultaneously
diagonalized if and only if the matrices commute. In the case of non-commuting
matrices, the best that can be achieved is simultaneous block diagonalization.
Here we give an efficient algorithm to explicitly compute a transfer matrix
which realizes the simultaneous block diagonalization of unitary matrices whose
decomposition in irreducible blocks (common invariant subspaces) is known from
elsewhere. Our main motivation lies in particle physics, where the resulting
transfer matrix must be known explicitly in order to unequivocally determine
the action of outer automorphisms such as parity, charge conjugation, or time
reversal on the particle spectrum.
