Fix MulStyle of CompositeMaps (#200)

Author: dkarrasch

src/composition.jl

@@ -19,6 +19,8 @@ Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::LinearMapTuple
 Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::AbstractVector{<:LinearMap{T}}) where {T} =
     CompositeMap{T}(reverse(maps))
 
+MulStyle(A::CompositeMap) = MulStyle(A.maps...) === TwoArg() ? TwoArg() : ThreeArg()
+
 # basic methods
 Base.size(A::CompositeMap) = (size(A.maps[end], 1), size(A.maps[1], 2))
 Base.axes(A::CompositeMap) = (axes(A.maps[end])[1], axes(A.maps[1])[2])
@@ -173,7 +175,7 @@ end
 
 function _unsafe_mul!(y, A::CompositeMap, x::AbstractVector)
     MulStyle(A) === TwoArg() ?
-        copyto!(y, foldr(*, reverse(A.maps), init=x)) :
+        copyto!(y, A*x) :
         _compositemul!(y, A, x)
     return y
 end

test/composition.jl

@@ -158,3 +158,28 @@ using LinearMaps: LinearMapVector, LinearMapTuple
         @test P * ones(3) == (LowerTriangular(ones(3,3))^i) * ones(3)
     end
 end
+
+# test product of 2-arg FunctionMaps
+# the following tests don't work when wrapped in a testset
+N = 100
+function planA()
+    y = zeros(N) # workspace
+    A = LinearMap{Float64}(x -> (y .= x .+ 1; y), N)
+    return A, y
+end
+function planB()
+    y = zeros(N) # workspace
+    A = LinearMap{Float64}(x -> (y .= x ./ 2; y), N)
+    return A, y
+end
+A, ya = planA()
+B, yb = planB()
+x = zeros(N)
+C = @inferred A*B; C*x
+@test C*x === ya == ones(N)
+D = @inferred B*A; D*x
+@test D*x === yb == fill(0.5, N)
+@test (@allocated C*x) == 0
+mul!(deepcopy(ya), C, x)
+y = deepcopy(ya)
+@test (@allocated mul!(y, C, x)) == 0