Fix type instability for liouvillian

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## Unreleased
 
 - Change the parameters structure of `sesolve`, `mesolve` and `mcsolve` functions to possibly support automatic differentiation. ([#311])
+- Fix type instability and reduce extra memory allocation in `liouvillian`. ([#315], [#318])
 
 
 ## [v0.21.5] (2024-11-15)
@@ -27,3 +28,5 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 [#306]: https://github.com/qutip/QuantumToolbox.jl/issues/306
 [#309]: https://github.com/qutip/QuantumToolbox.jl/issues/309
 [#311]: https://github.com/qutip/QuantumToolbox.jl/issues/311
+[#315]: https://github.com/qutip/QuantumToolbox.jl/issues/315
+[#318]: https://github.com/qutip/QuantumToolbox.jl/issues/318

src/qobj/superoperators.jl

@@ -24,15 +24,15 @@ end
 ## the rest of the SciMLOperators will just use lazy tensor (and prompt a warning)
 _spre(A::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spre(A.A, Id))
 _spre(A::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(A.λ, _spre(A.L, Id))
-_spre(A::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _spre(op, Id), +, A.ops)
+_spre(A::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _spre(op, Id), A.ops))
 function _spre(A::AbstractSciMLOperator, Id::AbstractMatrix)
     _lazy_tensor_warning("spre", A)
     return kron(Id, A)
 end
 
 _spost(B::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spost(B.A, Id))
 _spost(B::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(B.λ, _spost(B.L, Id))
-_spost(B::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _spost(op, Id), +, B.ops)
+_spost(B::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _spost(op, Id), B.ops))
 function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
     _lazy_tensor_warning("spost", B)
     return kron(transpose(B), Id)
@@ -43,7 +43,7 @@ _liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,Abstract
     -1im * (_spre(H, Id) - _spost(H, Id))
 _liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
 _liouvillian(H::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(H.λ, _liouvillian(H.L, Id))
-_liouvillian(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian(op, Id), +, H.ops)
+_liouvillian(H::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _liouvillian(op, Id), H.ops))
 
 # intrinsic lindblad_dissipator
 function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
@@ -166,7 +166,7 @@ function liouvillian(
 ) where {DT,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     L = liouvillian(H, Id_cache)
     if !(c_ops isa Nothing)
-        # sum all the Qobj first
+        # sum all the (time-independent) c_ops first
         c_ops_ti = filter(op -> isa(op, QuantumObject), c_ops)
         if !isempty(c_ops_ti)
             L += mapreduce(op -> lindblad_dissipator(op, Id_cache), +, c_ops_ti)

test/core-test/quantum_objects_evo.jl

@@ -213,7 +213,15 @@
         @test_throws ArgumentError cache_operator(L_td, ψ)
 
         @testset "Type Inference" begin
-            @inferred liouvillian(H_td, (a, QobjEvo(a', coef1)))
+            # we use destroy and create here because they somehow causes type instability before
+            H_td2 = H_td + QobjEvo(destroy(N) + create(N), coef3)
+            c_ops1 = (destroy(N), create(N))
+            c_ops2 = (destroy(N), QobjEvo(create(N), coef1))
+
+            @inferred liouvillian(H_td, c_ops1)
+            @inferred liouvillian(H_td, c_ops2)
+            @inferred liouvillian(H_td2, c_ops1)
+            @inferred liouvillian(H_td2, c_ops2)
         end
     end
 end