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