Computation of a simultaneous block-diagonalization

savecancel

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 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})$$ I know how 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