Implement generic operator type (#265)

* Switch from gem to mul

* Implement generalized Operator type

* Allow arbitrary array types as data

* Add some tests for sparse density matrices

* Loosen type constraints in dmaster_nh to work with lazy adjoints

* Add tests for irreducibility checks

* Better typing for alpha and beta in dmaster

* Infer types for alpha and beta in derivative functions

* Use eltype for types of alpha and beta in derivatives

* Better typing in stochastics

* Loosen type constraints on tspan

* Update version and dependencies

Project.toml

@@ -1,6 +1,6 @@
 name = "QuantumOptics"
 uuid = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
-version = "v0.7.3"
+version = "v0.8.0"
 
 [deps]
 QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
@@ -17,7 +17,7 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 
 [compat]
-QuantumOpticsBase = "0.1"
+QuantumOpticsBase = "0.2"
 Arpack = "0.4"
 DiffEqCallbacks = "2"
 FFTW = "1"
@@ -26,7 +26,7 @@ RecursiveArrayTools = "2"
 Reexport = "0.2"
 StochasticDiffEq = "6"
 WignerSymbols = "1"
-julia = "1"
+julia = "1.3"
 
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/bloch_redfield_master.jl

@@ -160,10 +160,10 @@ Time-evolution according to a Bloch-Redfield master equation.
         be changed.
 * `kwargs...`: Further arguments are passed on to the ode solver.
 """
-function master_bloch_redfield(tspan::Vector{Float64},
+function master_bloch_redfield(tspan,
         rho0::T, L::SuperOperator{Tuple{B,B},Tuple{B,B}},
         H::AbstractOperator{B,B}; fout::Union{Function,Nothing}=nothing,
-        kwargs...) where {B<:Basis,T<:DenseOperator{B,B}}
+        kwargs...) where {B<:Basis,T<:Operator{B,B}}
 
     #Prep basis transf
     evals, transf_mat = eigen(dense(H).data)
@@ -173,56 +173,58 @@ function master_bloch_redfield(tspan::Vector{Float64},
     # rho as Ket and L as DataOperator
     basis_comp = rho0.basis_l^2
     rho0_eb = Ket(basis_comp, (inv_transf_op * rho0 * transf_op).data[:]) #Transform to H eb and convert to Ket
-    L_ = isa(L, DenseSuperOperator) ? DenseOperator(basis_comp, L.data) : SparseOperator(basis_comp, L.data)
+    L_ = isa(L, SparseSuperOpType) ? SparseOperator(basis_comp, L.data) : DenseOperator(basis_comp, L.data)
 
     # Derivative function
-    dmaster_br_(t::Float64, rho::T2, drho::T2) where T2<:Ket = dmaster_br(drho, rho, L_)
+    dmaster_br_(t, rho::T2, drho::T2) where T2<:Ket = dmaster_br(drho, rho, L_)
 
     return integrate_br(tspan, dmaster_br_, rho0_eb, transf_op, inv_transf_op, fout; kwargs...)
 end
-master_bloch_redfield(tspan::Vector{Float64}, psi::Ket, args...; kwargs...) = master_bloch_redfield(tspan::Vector{Float64}, dm(psi), args...; kwargs...)
+master_bloch_redfield(tspan, psi::Ket, args...; kwargs...) = master_bloch_redfield(tspan, dm(psi), args...; kwargs...)
 
 # Derivative ∂ₜρ = Lρ
 function dmaster_br(drho::T, rho::T, L::DataOperator{B,B}) where {B<:Basis,T<:Ket{B}}
-    QuantumOpticsBase.gemv!(1.0, L, rho, 0.0, drho)
+    QuantumOpticsBase.mul!(drho,L,rho)
 end
 
 # Integrate if there is no fout specified
-function integrate_br(tspan::Vector{Float64}, dmaster_br::Function, rho::T,
+function integrate_br(tspan, dmaster_br::Function, rho::T,
                 transf_op::T2, inv_transf_op::T2, ::Nothing;
-                kwargs...) where {T<:Ket,T2<:DenseOperator}
+                kwargs...) where {T<:Ket,T2<:Operator}
     # Pre-allocate for in-place back-transformation from eigenbasis
     rho_out = copy(transf_op)
     tmp = copy(transf_op)
     tmp2 = copy(transf_op)
 
     # Define fout
-    function fout(t::Float64, rho::T)
+    function fout(t, rho::T)
         tmp.data[:] = rho.data
-        QuantumOpticsBase.gemm!(1.0, transf_op, tmp, 0.0, tmp2)
-        QuantumOpticsBase.gemm!(1.0, tmp2, inv_transf_op, 0.0, rho_out)
+        QuantumOpticsBase.mul!(tmp2,transf_op,tmp)
+        QuantumOpticsBase.mul!(rho_out,tmp2,inv_transf_op)
         return copy(rho_out)
     end
 
     return integrate(tspan, dmaster_br, copy(rho.data), rho, copy(rho), fout; kwargs...)
 end
 
 # Integrate with given fout
-function integrate_br(tspan::Vector{Float64}, dmaster_br::Function, rho::T,
+function integrate_br(tspan, dmaster_br::Function, rho::T,
                 transf_op::T2, inv_transf_op::T2, fout::Function;
-                kwargs...) where {T<:Ket,T2<:DenseOperator}
+                kwargs...) where {T<:Ket,T2<:Operator}
     # Pre-allocate for in-place back-transformation from eigenbasis
     rho_out = copy(transf_op)
     tmp = copy(transf_op)
     tmp2 = copy(transf_op)
 
+    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
+
     # Perform back-transfomration before calling fout
-    function fout_(t::Float64, rho::T)
+    function fout_(t, rho::T)
         tmp.data[:] = rho.data
-        QuantumOpticsBase.gemm!(1.0, transf_op, tmp, 0.0, tmp2)
-        QuantumOpticsBase.gemm!(1.0, tmp2, inv_transf_op, 0.0, rho_out)
+        QuantumOpticsBase.mul!(tmp2,transf_op,tmp)
+        QuantumOpticsBase.mul!(rho_out,tmp2,inv_transf_op)
         return fout(t, rho_out)
     end
 
-    return integrate(tspan, dmaster_br, copy(rho.data), rho, copy(rho), fout_; kwargs...)
+    return integrate(tspan_, dmaster_br, copy(rho.data), rho, copy(rho), fout_; kwargs...)
 end