Add herm_matmul to speed up mesolve.

[Ericgig]

Description
Add a variant to matmul specialised for case where the right matrix is a column stacked hermitian matrix, and the output is the same. Such as for mesolve using super operators.

Open solvers mesolve, brmesolve, smesolve and fmesolve all have a new options use_herm_matmul to use this operation. The default is False since we cannot easily check for exception (we can't test that a time dependant Hamiltonian is Hermitian, ...).

With large enough systems, this result in a visible speed up:

>>> N = 100
>>> H = qt.rand_herm(N, density=3/N)
>>> c_ops = [qt.destroy(N)]

>>> qt.mesolve(H, qt.basis(N, N-2), [0, 10], c_ops=c_ops, options={"use_herm_matmul": True}).stats["run time"]
0.8431086540222168
>>> qt.mesolve(H, qt.basis(N, N-2), [0, 10], c_ops=c_ops, options={"use_herm_matmul": False}).stats["run time"]
1.0559077262878418

qt.smesolve(H, qt.basis(N, N-2), [0, 2], sc_ops=c_ops, ntraj=1, options={"use_herm_matmul": True}).stats["run time"]
0.9043550491333008
qt.smesolve(H, qt.basis(N, N-2), [0, 2], sc_ops=c_ops, ntraj=1, options={"use_herm_matmul": False}).stats["run time"]
2.111673355102539

This is a use case for the capacity to dispatch on Data added in #2157. When herm_matmul is not available, it can be better to fallback on matmul using the same type than to do conversions between data types. For cupy, moving the data to the cpu to halves the work is probably not worth it.

Related issues or PRs
This is build on top of #2157, it should be merged first.

[coveralls]

coverage: 84.644% (-0.01%) from 84.655% when pulling 534cd4c on Ericgig:feature.herm_matmul into 5fa0ca6 on qutip:master.

[hodgestar]

I'd like to take a more detailed look at this once the Data specialization PR is merged, but the concept looks good.

Would you mind pinging me for another review here once the other PR is done?

[hodgestar]

I left a couple more comments. On the whole I'm not wild about the extra QobjEvo class, but if this is a use case you'd like to support, what about the following suggestion:

We allow the user to supply the QobjEvo class to use as an option, and allow them to specify either the class itself or a name for the class. E.g. qobjevo_cls=QobjEvoHerm or qobjevo_cls=herm.

That allows many different classes to be used without the solvers needing to have special support added each time.

Perhaps we should also not allow users to switch the class later by changing options? That might simplify the logic because we wouldn't need to modify the RHS when options are updated.

[hodgestar]

Are these two lines meant to be duplicated here? It seems odd to call self._update_rhs() twice.

[hodgestar]

Duplicating the RHS here feels very messy. Could we just store a single RHS as self.rhs and convert to and from QobjEvoHerm as needed when options are updated?

[Ericgig]

We allow the user to supply the QobjEvo class to use as an option, and allow them to specify either the class itself or a name for the class. E.g. qobjevo_cls=QobjEvoHerm or qobjevo_cls=herm.

I would prefer the user not needing to know how we do it, just that there an option that speed up the simulation by 40 % in normal cases. Knowing that we forced it in an alternative qobjevo class is not useful and I hope we won't be forced to add many kinds of qobjevo.

Perhaps we should also not allow users to switch the class later by changing options? That might simplify the logic because we wouldn't need to modify the RHS when options are updated.

Since changing options does not change the physic, I would like them to be changeable. But the rhs, and _rhs is certainly not great. I will rethink the way to do it.

[Ericgig]

@hodgestar

I removed the _rhs by changing it so the rhs is built when making the integrator instead of in __init__.

I also simplified the way options are updated so brmesolve no longer need to overwrite _apply_options.

[hodgestar]

All subclasses of Solver should have one path for building the RHS. Having two paths makes the code much harder to reason about.

I think _build_rhs() is a good idea, although I think we should implement it on all sub-classes of Solver and have Solver do self.rhs = self._build_rhs() at the end of __init__ and when options are changed (if needed). That way .rhs is always correct for the given options.