Store domainspace in TimesOperator

Author: jishnub

src/Operators/general/algebra.jl

@@ -223,19 +223,23 @@ axpy!(α, S::SubOperator{T,OP}, A::AbstractMatrix) where {T,OP<:ConstantTimesOpe
 
 
 
-struct TimesOperator{T,BW,SZ,O<:Operator{T},BBW,SBBW} <: Operator{T}
+struct TimesOperator{T,BW,SZ,O<:Operator{T},BBW,SBBW,D} <: Operator{T}
     ops::Vector{O}
     bandwidths::BW
     sz::SZ
     blockbandwidths::BBW
     subblockbandwidths::SBBW
     isbandedblockbanded::Bool
     israggedbelow::Bool
+    domainspace::D
 
     function TimesOperator{T,BW,SZ,O,BBW,SBBW}(ops::Vector{O}, bw::BW,
         sz::SZ, bbw::BBW, sbbw::SBBW,
-        ibbb::Bool, irb::Bool) where {T,O<:Operator{T},BW,SZ,BBW,SBBW}
-        # check compatible
+        ibbb::Bool, irb::Bool, dsp::D) where {T,O<:Operator{T},BW,SZ,BBW,SBBW,D}
+
+        dsp == domainspace(ops[end]) || throw(ArgumentError("incompatible domainspace"))
+
+        # check compatibility
         for k = 1:length(ops)-1
             size(ops[k], 2) == size(ops[k+1], 1) || throw(ArgumentError("incompatible operator sizes"))
             spacescompatible(domainspace(ops[k]), rangespace(ops[k+1])) || throw(ArgumentError("incompatible spaces at index $k"))
@@ -252,7 +256,7 @@ struct TimesOperator{T,BW,SZ,O<:Operator{T},BBW,SBBW} <: Operator{T}
             newops = ops
         end
 
-        new{T,BW,SZ,O,BBW,SBBW}(newops, bw, sz, bbw, sbbw, ibbb, irb)
+        new{T,BW,SZ,O,BBW,SBBW,D}(newops, bw, sz, bbw, sbbw, ibbb, irb, dsp)
     end
 end
 
@@ -273,9 +277,10 @@ function TimesOperator(ops::AbstractVector{O},
         sbbw::Tuple{Any,Any}=bandwidthssum(ops, subblockbandwidths),
         ibbb::Bool=all(isbandedblockbanded, ops),
         irb::Bool=all(israggedbelow, ops),
+        dsp = domainspace(last(ops)),
         ) where {O<:Operator}
     TimesOperator{eltype(O),typeof(bw),typeof(sz),O,typeof(bbw),typeof(sbbw)}(
-        convert_vector(ops), bw, sz, bbw, sbbw, ibbb, irb)
+        convert_vector(ops), bw, sz, bbw, sbbw, ibbb, irb, dsp)
 end
 
 _extractops(A::TimesOperator, ::typeof(*)) = A.ops
@@ -284,9 +289,10 @@ function TimesOperator(A::Operator, B::Operator)
     v = collateops(*, A, B)
     ibbb = all(isbandedblockbanded, (A, B))
     irb = all(israggedbelow, (A, B))
+    dsp = domainspace(B)
     TimesOperator(convert_vector(v), _bandwidthssum(A, B), _timessize((A, B)),
         _bandwidthssum(A, B, blockbandwidths),
-        _bandwidthssum(A, B, subblockbandwidths), ibbb, irb)
+        _bandwidthssum(A, B, subblockbandwidths), ibbb, irb, dsp)
 end
 
 
@@ -301,7 +307,7 @@ function convert(::Type{Operator{T}}, P::TimesOperator) where {T}
                       _convertops(Operator{T}, ops),
             bandwidths(P), size(P), blockbandwidths(P),
             subblockbandwidths(P), isbandedblockbanded(P),
-            israggedbelow(P))::Operator{T}
+            israggedbelow(P), domainspace(P))::Operator{T}
     end
 end
 
@@ -318,7 +324,7 @@ end
     @assert length(opsin) > 1 "need at least 2 operators"
     ops, bw, bbw, sbbw, ibbb, irb = __promotetimes(opsin, dsp, anytimesop)
     sz = _timessize(ops)
-    TimesOperator(convert_vector(ops), bw, sz, bbw, sbbw, ibbb, irb)
+    TimesOperator(convert_vector(ops), bw, sz, bbw, sbbw, ibbb, irb, dsp)
 end
 function __promotetimes(opsin, dsp, anytimesop)
     ops = Vector{Operator{promote_eltypeof(opsin)}}(undef, 0)
@@ -372,7 +378,8 @@ end
     end
 end
 
-domainspace(P::PlusOrTimesOp) = domainspace(last(P.ops))
+domainspace(P::PlusOperator) = domainspace(last(P.ops))
+domainspace(T::TimesOperator) = T.domainspace
 rangespace(P::PlusOrTimesOp) = rangespace(first(P.ops))
 
 domain(P::TimesOperator) = commondomain(P.ops)

test/runtests.jl

@@ -259,9 +259,10 @@ end
             f = Fun(sp, coeff)
             for sp2 in Any[(), (sp,)]
                 M = Multiplication(f, sp2...)
-                a = (M * M) * M
-                b = M * (M * M)
+                a = TimesOperator(M, M) * M
+                b = M * TimesOperator(M, M)
                 @test a == b
+                @test (@inferred domainspace(a)) == domainspace(M)
                 @test bandwidths(a) == bandwidths(b)
             end
             M = Multiplication(Fun(PointSpace(1:3)))