From 19a9d53f115991fd9e29204012bf7de3b1e91f31 Mon Sep 17 00:00:00 2001
From: Akira Kyle <akira@akirakyle.com>
Date: Tue, 12 Sep 2023 14:38:40 -0600
Subject: [PATCH] Fix sparse matrix exponential of zero since FastExpm doesn't
 handle it (#136)

---
 src/operators_sparse.jl     | 6 +++++-
 src/superoperators.jl       | 9 ++++++++-
 test/test_operators.jl      | 1 +
 test/test_superoperators.jl | 1 +
 4 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/operators_sparse.jl b/src/operators_sparse.jl
index 7c09dd75..b52d5ce7 100644
--- a/src/operators_sparse.jl
+++ b/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}
diff --git a/src/superoperators.jl b/src/superoperators.jl
index 3a137752..8128c3e5 100644
--- a/src/superoperators.jl
+++ b/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
diff --git a/test/test_operators.jl b/test/test_operators.jl
index 92fc8436..5e180441 100644
--- a/test/test_operators.jl
+++ b/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]))
diff --git a/test/test_superoperators.jl b/test/test_superoperators.jl
index 28043bdb..204264ac 100644
--- a/test/test_superoperators.jl
+++ b/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