Fix sparse matrix exponential of zero since FastExpm doesn't handle it (

#136)

src/operators_sparse.jl

@@ -56,7 +56,11 @@ If you only need the result of the exponential acting on a vector,
 consider using much faster implicit methods that do not calculate the entire exponential.
 """
 function exp(op::T; opts...) where {B,T<:SparseOpType{B,B}}
-    return SparseOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...))
+    if iszero(op)
+        return identityoperator(op)
+    else
+        return SparseOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...))
+    end
 end
 
 function permutesystems(rho::SparseOpPureType{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis}

src/superoperators.jl

@@ -241,9 +241,16 @@ All optional arguments are passed to `fastExpm` and can be used to specify toler
 If you only need the result of the exponential acting on an operator,
 consider using much faster implicit methods that do not calculate the entire exponential.
 """
-Base.exp(op::SparseSuperOpType; opts...) = SuperOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...))
+function Base.exp(op::SparseSuperOpType; opts...)
+    if iszero(op)
+        return identitysuperoperator(op)
+    else
+        return SuperOperator(op.basis_l, op.basis_r, fastExpm(op.data; opts...))
+    end
+end
 
 # Array-like functions
+Base.zero(A::SuperOperator) = SuperOperator(A.basis_l, A.basis_r, zero(A.data))
 Base.size(A::SuperOperator) = size(A.data)
 @inline Base.axes(A::SuperOperator) = axes(A.data)
 Base.ndims(A::SuperOperator) = 2

test/test_operators.jl

@@ -118,6 +118,7 @@ b = FockBasis(40)
 alpha = 1+5im
 H = alpha * create(b) - conj(alpha) * destroy(b)
 @test exp(sparse(H); threshold=1e-10) ≈ displace(b, alpha)
+@test exp(sparse(zero(identityoperator(b)))) ≈ identityoperator(b)
 
 @test one(b1).data == Diagonal(ones(b1.shape[1]))
 @test one(op1).data == Diagonal(ones(b1.shape[1]))

test/test_superoperators.jl

@@ -221,6 +221,7 @@ Ldense .+= L
 b = FockBasis(20)
 L = liouvillian(identityoperator(b), [destroy(b)])
 @test exp(sparse(L)).data ≈ exp(dense(L)).data
+@test exp(sparse(zero(identitysuperoperator(b)))).data ≈ identitysuperoperator(b).data
 N = exp(log(2) * sparse(L)) # 50% loss channel
 ρ = N * dm(coherentstate(b, 1))
 @test (1 - real(tr(ρ^2))) < 1e-10 # coherent state remains pure under loss