Mul adjoint and lazy sum (#85)

* allocate concrete type operator in *

* tests

* fix test

* Robust version of _parent

* Add test with sparse sparse mul

---------

Co-authored-by: Amit Rotem <[email protected]>

Author: AmitRotem

src/operators_dense.jl

@@ -89,12 +89,12 @@ Base.isapprox(x::DataOperator, y::DataOperator; kwargs...) = false
 *(a::Operator, b::Number) = Operator(a.basis_l, a.basis_r, b*a.data)
 *(a::Number, b::Operator) = Operator(b.basis_l, b.basis_r, a*b.data)
 function *(op1::AbstractOperator{B1,B2}, op2::Operator{B2,B3,T}) where {B1,B2,B3,T}
-    result = Operator{B1,B3,T}(op1.basis_l, op2.basis_r, similar(op2.data,length(op1.basis_l),length(op2.basis_r)))
+    result = Operator{B1,B3}(op1.basis_l, op2.basis_r, similar(_parent(op2.data),promote_type(eltype(op1),eltype(op2)),length(op1.basis_l),length(op2.basis_r)))
     mul!(result,op1,op2)
     return result
 end
 function *(op1::Operator{B1,B2,T}, op2::AbstractOperator{B2,B3}) where {B1,B2,B3,T}
-    result = Operator{B1,B3,T}(op1.basis_l, op2.basis_r, similar(op1.data,length(op1.basis_l),length(op2.basis_r)))
+    result = Operator{B1,B3}(op1.basis_l, op2.basis_r, similar(_parent(op1.data),promote_type(eltype(op1),eltype(op2)),length(op1.basis_l),length(op2.basis_r)))
     mul!(result,op1,op2)
     return result
 end
@@ -109,6 +109,10 @@ function *(psi::Bra{BL,T}, op::AbstractOperator{BL,BR}) where {BL,BR,T}
     return result
 end
 
+_parent(x::T, x_parent::T) where T = x
+_parent(x, x_parent) = _parent(x_parent, parent(x_parent))
+_parent(x) = _parent(x, parent(x))
+
 /(a::Operator, b::Number) = Operator(a.basis_l, a.basis_r, a.data ./ b)
 
 dagger(x::Operator) = Operator(x.basis_r,x.basis_l,adjoint(x.data))

test/test_operators_lazysum.jl

@@ -102,6 +102,25 @@ 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_))
 
+## multiplication with Operator of AbstractMatrix
+LSop = LazySum(randoperator(b1a^2)) # AbstractOperator
+LSop_s = LazySum(sparse(randoperator(b1a^2)))
+hermitian_op = Operator(basis(LSop), Hermitian(randn(ComplexF64,length(basis(LSop)),length(basis(LSop))))) # Hermitian
+symmetric_op = Operator(basis(LSop), Symmetric(randn(ComplexF64,length(basis(LSop)),length(basis(LSop))))) # Symmetric
+adjoint_op = randoperator(basis(LSop))' # Adjoint
+real_op = Operator(basis(LSop), real(adjoint_op.data)) # real
+ops_tuple = (symmetric_op,hermitian_op,adjoint_op)
+### Test
+@test ops_tuple.*(LSop,) == dense.(ops_tuple).*(LSop,)
+@test (LSop,).*ops_tuple == (LSop,).*dense.(ops_tuple)
+### test with sparse
+@test all(isapprox.(sparse.(ops_tuple).*(LSop,) , dense.(ops_tuple).*(LSop,), atol=1e-13))
+@test all(isapprox.((LSop,).*sparse.(ops_tuple) , (LSop,).*dense.(ops_tuple), atol=1e-13))
+@test all(isapprox.(dense.(sparse.(ops_tuple).*(LSop_s,)) , dense.(ops_tuple).*(LSop_s,), atol=1e-13))
+@test all(isapprox.(dense.((LSop_s,).*sparse.(ops_tuple)) , (LSop_s,).*dense.(ops_tuple), atol=1e-13))
+### test real valued op with AbstractOperator
+@test isapprox(LSop*real_op*LSop , LSop*Operator(basis(real_op), complex.(real_op.data))*LSop, atol=1e-13)
+
 # Test division
 @test 1e-14 > D(op1/7, op1_/7)