Fix for multiplication of LazySum with no elements (#95)

* eltype LazySum with no elements

* test multiplication with LazySum with no elements

* add zero LazySum multiplication test

* add dim check for LazySum with no elements

* more tests for zero LazySum

* small fix

* zero method for Bra and Ket

Author: AmitRotem

src/operators_lazysum.jl

@@ -59,7 +59,7 @@ LazySum(operators::AbstractOperator...) = LazySum(mapreduce(eltype, promote_type
 LazySum() = throw(ArgumentError("LazySum needs a basis, or at least one operator!"))
 
 Base.copy(x::LazySum) = @samebases LazySum(x.basis_l, x.basis_r, copy(x.factors), copy.(x.operators))
-Base.eltype(x::LazySum) = promote_type(eltype(x.factors), mapreduce(eltype, promote_type, x.operators))
+Base.eltype(x::LazySum) = mapreduce(eltype, promote_type, x.operators; init=eltype(x.factors))
 
 dense(op::LazySum) = length(op.operators) > 0 ? sum(op.factors .* dense.(op.operators)) : Operator(op.basis_l, op.basis_r, zeros(eltype(op.factors), length(op.basis_l), length(op.basis_r)))
 SparseArrays.sparse(op::LazySum) = length(op.operators) > 0 ? sum(op.factors .* sparse.(op.operators)) : Operator(op.basis_l, op.basis_r, spzeros(eltype(op.factors), length(op.basis_l), length(op.basis_r)))
@@ -140,6 +140,7 @@ end
 # Fast in-place multiplication
 function mul!(result::Ket{B1},a::LazySum{B1,B2},b::Ket{B2},alpha,beta) where {B1,B2}
     if length(a.operators) == 0
+        _check_mul!_dim_compatibility(size(result), size(a), size(b))
         result.data .*= beta
     else
         mul!(result,a.operators[1],b,alpha*a.factors[1],beta)
@@ -152,6 +153,7 @@ end
 
 function mul!(result::Bra{B2},a::Bra{B1},b::LazySum{B1,B2},alpha,beta) where {B1,B2}
     if length(b.operators) == 0
+        _check_mul!_dim_compatibility(size(result), reverse(size(b)), size(a))
         result.data .*= beta
     else
         mul!(result,a,b.operators[1],alpha*b.factors[1],beta)
@@ -164,6 +166,7 @@ end
 
 function mul!(result::Operator{B1,B3},a::LazySum{B1,B2},b::Operator{B2,B3},alpha,beta) where {B1,B2,B3}
     if length(a.operators) == 0
+        _check_mul!_dim_compatibility(size(result), size(a), size(b))
         result.data .*= beta
     else
         mul!(result,a.operators[1],b,alpha*a.factors[1],beta)
@@ -175,6 +178,7 @@ function mul!(result::Operator{B1,B3},a::LazySum{B1,B2},b::Operator{B2,B3},alpha
 end
 function mul!(result::Operator{B1,B3},a::Operator{B1,B2},b::LazySum{B2,B3},alpha,beta) where {B1,B2,B3}
     if length(b.operators) == 0
+        _check_mul!_dim_compatibility(size(result), size(a), size(b))
         result.data .*= beta
     else
         mul!(result,a,b.operators[1],alpha*b.factors[1],beta)

src/states.jl

@@ -30,6 +30,8 @@ mutable struct Ket{B,T} <: AbstractKet{B,T}
     end
 end
 
+Base.zero(x::Bra) = Bra(x.basis, zero(x.data))
+Base.zero(x::Ket) = Ket(x.basis, zero(x.data))
 eltype(::Type{K}) where {K <: Ket{B,V}} where {B,V} = eltype(V)
 eltype(::Type{K}) where {K <: Bra{B,V}} where {B,V} = eltype(V)
 

test/test_operators_lazysum.jl

@@ -102,6 +102,16 @@ xbra1 = Bra(b_l, rand(ComplexF64, length(b_l)))
 @test 1e-11 > D((op1+op2)*(x1+0.3*x2), (op1_+op2_)*(x1+0.3*x2))
 @test 1e-12 > D(dagger(x1)*dagger(0.3*op2), dagger(x1)*dagger(0.3*op2_))
 
+## Test multiplication with LazySum that has no elements
+@test iszero( LazySum(b_r, b_l) * op1a )
+@test iszero( op1a * LazySum(b_r, b_l) )
+@test iszero( LazySum(b_l, b_r) * x1 )
+@test iszero( xbra1 * LazySum(b_l, b_r) )
+@test_throws DimensionMismatch LazySum(FockBasis(2), NLevelBasis(2)) * randoperator(NLevelBasis(4), GenericBasis(2)) # save Basis with different size
+@test_throws DimensionMismatch randoperator(GenericBasis(1), FockBasis(3)) * LazySum(FockBasis(1), NLevelBasis(2))
+@test_throws DimensionMismatch LazySum(FockBasis(2), NLevelBasis(2)) * randstate(NLevelBasis(7))
+@test_throws DimensionMismatch randstate(FockBasis(3))' * LazySum(FockBasis(1), NLevelBasis(2))
+
 ## multiplication with Operator of AbstractMatrix
 LSop = LazySum(randoperator(b1a^2)) # AbstractOperator
 LSop_s = LazySum(sparse(randoperator(b1a^2)))