Fix bug in TimesOperator constructor, and add Multiplication equality (#110)

* performance optimizations in TimesOperator

* add equality methods and tests

Author: jishnub

src/Operators/banded/Multiplication.jl

@@ -26,6 +26,8 @@ function ConcreteMultiplication(f::Fun{D,T},sp::Space) where {D,T}
     ConcreteMultiplication{D,typeof(sp),V}(convert(Fun{D,V},chop(f,40*eps(cfstype(f)))),sp)
 end
 
+==(A::ConcreteMultiplication, B::ConcreteMultiplication) = (A.f == B.f) && (A.space == B.space)
+
 # We do this in two stages to support Modifier spaces
 # without ambiguity errors
 function defaultMultiplication(f::Fun,sp::Space)

src/Operators/general/FiniteOperator.jl

@@ -19,6 +19,7 @@ FiniteOperator(M::AbstractMatrix{<:Number}) =
 convert(::Type{Operator{T}},F::FiniteOperator) where {T} =
     FiniteOperator(convert(AbstractMatrix{T},F.matrix),F.domainspace,F.rangespace)::Operator{T}
 
+==(A::FiniteOperator, B::FiniteOperator) = A.matrix == B.matrix && A.domainspace == B.domainspace && A.rangespace == B.rangespace
 
 Base.promote_rule(::Type{OT},::Type{MT}) where {OT<:Operator,MT<:AbstractMatrix} = Operator{promote_type(eltype(OT),eltype(MT))}
 

src/Operators/general/algebra.jl

@@ -228,45 +228,27 @@ struct TimesOperator{T,BI} <: Operator{T}
     function TimesOperator{T,BI}(ops::Vector{Operator{T}},bi::BI) where {T,BI}
         # check compatible
         for k=1:length(ops)-1
-            @assert size(ops[k],2) == size(ops[k+1],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("imcompatible spaces at index $k"))
         end
 
         # remove TimesOperators buried inside ops
-        hastimes = false
-        for k=1:length(ops)-1
-            @assert spacescompatible(domainspace(ops[k]),rangespace(ops[k+1]))
-            hastimes = hastimes || isa(ops[k],TimesOperator)
-        end
-
-        if hastimes
-            newops=Vector{Operator{T}}(0)
-            for op in ops
-               if isa(op,TimesOperator)
-                    for op2 in op.ops
-                        push!(newops,op2)
-                    end
-                else
-                    push!(newops,op)
-                end
+        timesinds = findall(x -> isa(x, TimesOperator), ops)
+        if !isempty(timesinds)
+            newops = copy(ops)
+            for ind in timesinds
+                splice!(newops, ind, ops[ind].ops)
             end
-            ops=newops
+        else
+            newops = ops
         end
 
-
-        new{T,BI}(ops,bi)
+        new{T,BI}(newops,bi)
     end
 end
 
-
-function bandwidthsum(P,k)
-    ret=0
-    for op in P
-        ret+=bandwidths(op)[k]
-    end
-    ret
-end
-
-bandwidthssum(P) = (bandwidthsum(P,1),bandwidthsum(P,2))
+bandwidthssum(P, k) = bandwidthssum(P)[k]
+bandwidthssum(P) = mapreduce(bandwidths, (t1, t2) -> t1 .+ t2, P, init = (0,0))
 
 TimesOperator(ops::Vector{Operator{T}},bi::Tuple{N1,N2}) where {T,N1,N2} =
     TimesOperator{T,typeof(bi)}(ops,bi)

test/runtests.jl

@@ -82,4 +82,36 @@ end
     @test ApproxFunBase.coefficients(f) === v
 end
 
+@testset "operator algebra" begin
+    @testset "Multiplication" begin
+        sp = PointSpace(1:3)
+        coeff = [1:3;]
+        f = Fun(sp, coeff)
+        for sp2 in Any[(), (sp,)]
+            a = Multiplication(f, sp2...)
+            b = Multiplication(f, sp2...)
+            @test a == b
+            @test bandwidths(a) == bandwidths(b)
+        end
+    end
+    @testset "TimesOperator" begin
+        sp = PointSpace(1:3)
+        coeff = [1:3;]
+        f = Fun(sp, coeff)
+        for sp2 in Any[(), (sp,)]
+            M = Multiplication(f, sp2...)
+            a = (M * M) * M
+            b = M * (M * M)
+            @test a == b
+            @test bandwidths(a) == bandwidths(b)
+        end
+        @testset "unwrap TimesOperator" begin
+            M = Multiplication(f)
+            for ops in Any[Operator{Float64}[M, M * M], Operator{Float64}[M*M, M]]
+                @test TimesOperator(ops).ops == [M, M, M]
+            end
+        end
+    end
+end
+
 @time include("ETDRK4Test.jl")