### Computation of a ~~simultenous ~~simultanenous block-diagonalization

Let $n$ be a positive integer and consider of finite set $S \subset M_n(\mathbb{C})$ such that $S^* = S$ (i.e. if $M \in S$ then $M^* \in S$). The algebra generated by $S$ is a finite dimensional $*$-algebra over $\mathbb{C}$, so is isomorphic to a direct sum of matrix ~~algebra (it is a particular case of Wedderburn decomposition), ~~algebra, 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})$$
I know I to get $\langle S \rangle$ by using `FiniteDimensionalAlgebra(CC,[M for M in S])`

.

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

*Remark*: When the matrices commute over each other, it is called a *simultaneous diagonalization*, and I know how to compute it using `jordan_form(transformation=True)`

several times.

In some sense, what I am looking for in general is how to compute a *simultaneous block-diagonalization*