Let $n$ be a positive integer and consider of finite set $S \subset M_n(\mathbb{C})$ such that $S^* = S$ (i.e. if $a \in S$ then $a^* \in S$). The algebra generated by $S$ is a finite dimensional $*$-algebra over $\mathbb{C}$, so is isomorphic to a direct sum of matrix algebras, i.e. there are $n_1 \le n_2 \le \dots \le n_r$ such that:

$$\langle S \rangle \simeq \bigoplus_{i=1}^r M_{n_i}(\mathbb{C})$$

**Question**: How to compute the change of basis $p$ such that for all $a \in S$, we have $p^{-1}ap$ block-diagonal as for the above decomposition?

We can see each $M_{n_i}(\mathbb{C})$ above as $End(V_i)$ with $V_i$ an irreducible representation of $\langle S \rangle$.

If the matrices commute over each other then (using $S^*=S$) we get that the matrices are normal, so diagonalizable, and then what we ask in this case is just a simultaneous diagonalization (ok).
So in general, what I am looking for is how to compute a (thinnest) simultaneous *block*-diagonalization.