Use diagonalization instead of scipy.expm in the cases where it's faster (matrix dim <~ 1600)

[miamico]

What is the expected behavior?

When taking the matrix exponential of an Hermitean matrix $M$ (for example to compute the time evolution operator from the Hamiltonian of a closed system) it is sometimes faster to first diagonalize the matrix $M = V D V^{dagger}$ and then exponentiate the diagonal matrix $D$ (which is given by exponentiating the diagonal elements) instead of using the built-in method in scipy expm directly on the matrix $M$.

Since fixed-step solvers are meant to work for arbitrary generators (even ones that cannot be diagonalized) but this method is faster for diagonalizable generators (in particular anti-hermitian), then it would be worth adding this implementation as an option for the user and to choose this option automatically in some higher level interfaces (e.g. when using a Solver instance that only has Hamiltonian information).

Possible implementation details:
Add the option to the core solver method, and add a check in the solve_lmde interface that if what's being simulated is a HamiltonianModel object, and the dim of the initial state is below some heuristically-set level, then it automatically uses the diagonalization method.