Merge pull request #2408 from Ericgig/solver.faster_normalizations

Faster dm norm in solvers

qutip/solver/solver_base.py

@@ -85,6 +85,9 @@ def _prepare_state(self, state):
 
         self._state_metadata = {
             'dims': state._dims,
+            # This is herm flag take for granted that the liouvillian keep
+            # hermiticity.  But we do not check user passed super operator for
+            # anything other than dimensions.
             'isherm': state.isherm and not (self.rhs.dims == state.dims)
         }
         if self.rhs.dims[1] == state.dims:
@@ -102,7 +105,10 @@ def _restore_state(self, data, *, copy=True):
             state = Qobj(data, **self._state_metadata, copy=copy)
 
         if data.shape[1] == 1 and self._options['normalize_output']:
-            state = state * (1 / state.norm())
+            if state.isoper:
+                state = state * (1 / state.tr())
+            else:
+                state = state * (1 / state.norm())
 
         return state
 

qutip/tests/solver/test_mesolve.py

@@ -206,7 +206,7 @@ def testME_TDDecayliouvillian(self, c_ops):
     def test_mesolve_normalization(self, state_type):
         # non-hermitean H causes state to evolve non-unitarily
         H = qutip.Qobj([[1, -0.1j], [-0.1j, 1]])
-        H = qutip.sprepost(H, H) # ensure use of MeSolve
+        H = qutip.spre(H) + qutip.spost(H.dag()) # ensure use of MeSolve
         psi0 = qutip.basis(2, 0)
         options = {"normalize_output": True, "progress_bar": None}