This does not answer your question and has nothing to do with scipy, but we have a somewhat specialized implementation of simultaneous diagonalization of hermitian matrices in qutip (a package for quantum dynamics)

Perhaps it could be useful for you when you work on your implementation. Also, if you have, or find, a better or more efficient way of implementing this function we (qutip project) would be interested in hearing from you as well.

Best regards

Robert