Slightly reorganize composition, improve coverage

Author: dkarrasch

src/composition.jl

@@ -22,8 +22,8 @@ Base.mapreduce(::typeof(identity), ::typeof(Base.mul_prod), maps::AbstractVector
 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])
+Base.size(A::CompositeMap) = (size(last(A.maps), 1), size(first(A.maps), 2))
+Base.axes(A::CompositeMap) = (axes(last(A.maps))[1], axes(first(A.maps))[2])
 Base.isreal(A::CompositeMap) = all(isreal, A.maps) # sufficient but not necessary
 
 # the following rules are sufficient but not necessary
@@ -32,17 +32,17 @@ for (f, _f, g) in ((:issymmetric, :_issymmetric, :transpose),
     @eval begin
         LinearAlgebra.$f(A::CompositeMap) = $_f(A.maps)
         $_f(maps::Tuple{}) = true
-        $_f(maps::Tuple{<:LinearMap}) = $f(maps[1])
+        $_f(maps::Tuple{<:LinearMap}) = $f(first(maps))
         $_f(maps::LinearMapTuple) =
-            maps[end] == $g(maps[1]) && $_f(Base.front(Base.tail(maps)))
+            maps[end] == $g(first(maps)) && $_f(Base.front(Base.tail(maps)))
         function $_f(maps::LinearMapVector)
             n = length(maps)
             if n == 0
                 return true
             elseif n == 1
-                return ($f(maps[1]))::Bool
+                return ($f(first(maps)))::Bool
             else
-                return ((maps[end] == $g(maps[1]))::Bool && $_f(@views maps[2:end-1]))
+                return ((last(maps) == $g(first(maps)))::Bool && $_f(@views maps[begin+1:end-1]))
             end
         end
         # since the introduction of ScaledMap, the following cases cannot occur
@@ -71,10 +71,10 @@ function _isposdef(maps::LinearMapVector)
     if n == 0
         return true
     elseif n == 1
-        return isposdef(maps[1])
+        return isposdef(first(maps))
     else
-        return (maps[end] == adjoint(maps[1]) || maps[end] == maps[1]) &&
-            isposdef(maps[1]) && _isposdef(maps[2:end-1])
+        return (last(maps) == adjoint(first(maps)) || last(maps) == first(maps)) &&
+            isposdef(first(maps)) && _isposdef(maps[begin+1:end-1])
     end
 end
 
@@ -180,12 +180,11 @@ function _unsafe_mul!(y, A::CompositeMap, x::AbstractVector)
     return y
 end
 _unsafe_mul!(y, A::CompositeMap, x::AbstractMatrix) = _compositemul!(y, A, x)
+_unsafe_mul!(y, A::CompositeMap{<:Any,<:Tuple{LinearMap}}, x::AbstractVector) =
+    _unsafe_mul!(y, A.maps[1], x)
+_unsafe_mul!(y, A::CompositeMap{<:Any,<:Tuple{LinearMap}}, X::AbstractMatrix) =
+    _unsafe_mul!(y, A.maps[1], X)
 
-function _compositemul!(y, A::CompositeMap{<:Any,<:Tuple{LinearMap}}, x,
-                        source = nothing,
-                        dest = nothing)
-    return _unsafe_mul!(y, A.maps[1], x)
-end
 function _compositemul!(y, A::CompositeMap{<:Any,<:Tuple{LinearMap,LinearMap}}, x,
                         source = nothing,
                         dest = nothing)
@@ -206,9 +205,9 @@ function _compositemul!(y, A::CompositeMap{<:Any,<:LinearMapVector}, x,
                         dest = nothing)
     N = length(A.maps)
     if N == 1
-        return _unsafe_mul!(y, A.maps[1], x)
+        return _unsafe_mul!(y, A.maps[begin], x)
     elseif N == 2
-        return _unsafe_mul!(y, A.maps[2] * A.maps[1], x)
+        return _unsafe_mul!(y, A.maps[end] * A.maps[begin], x)
     else
         return _compositemulN!(y, A, x, source, dest)
     end
@@ -218,19 +217,21 @@ function _compositemulN!(y, A::CompositeMap, x,
                          src = nothing,
                          dst = nothing)
     N = length(A.maps) # ≥ 3
+    n = n0 = firstindex(A.maps)
     source = isnothing(src) ?
-        convert(AbstractArray, A.maps[1] * x) :
-        _unsafe_mul!(src, A.maps[1], x)
+        convert(AbstractArray, A.maps[n] * x) :
+        _unsafe_mul!(src, A.maps[n], x)
+    n += 1
     dest = isnothing(dst) ?
-        convert(AbstractArray, A.maps[2] * source) :
-        _unsafe_mul!(dst, A.maps[2], source)
+        convert(AbstractArray, A.maps[n] * source) :
+        _unsafe_mul!(dst, A.maps[n], source)
     dest, source = source, dest # alternate dest and source
-    for n in 3:N-1
+    for n in (n0+2):N-1
         dest = _resize(dest, (size(A.maps[n], 1), size(x)[2:end]...))
         _unsafe_mul!(dest, A.maps[n], source)
         dest, source = source, dest # alternate dest and source
     end
-    _unsafe_mul!(y, A.maps[N], source)
+    _unsafe_mul!(y, last(A.maps), source)
     return y
 end
 

test/composition.jl

@@ -6,6 +6,10 @@ using LinearMaps: LinearMapVector, LinearMapTuple
     @test F == LinearMap(cumsum, reverse ∘ cumsum ∘ reverse, 10; ismutating=false)
     FC = LinearMap{ComplexF64}(cumsum, reverse ∘ cumsum ∘ reverse, 10; ismutating=false)
     FCM = @inferred LinearMaps.CompositeMap{ComplexF64}((FC,))
+    FCMv = @inferred LinearMaps.CompositeMap{ComplexF64}([FC,])
+    FCiip = LinearMaps.CompositeMap{ComplexF64}((LinearMap{ComplexF64}(cumsum!, 10),))
+    FCM2v = @inferred LinearMaps.CompositeMap{ComplexF64}([FC, FC])
+    FCM2iipv = @inferred LinearMaps.CompositeMap{ComplexF64}([FCiip, FCiip])
     L = LowerTriangular(ones(10,10))
     @test_throws DimensionMismatch F * LinearMap(zeros(2,2))
     @test_throws ErrorException LinearMaps.CompositeMap{Float64}((FC, LinearMap(rand(10,10))))
@@ -17,7 +21,8 @@ using LinearMaps: LinearMapVector, LinearMapTuple
     N = @inferred LinearMap(B)
     v = rand(ComplexF64, 10)
     α = rand(ComplexF64)
-    @test FCM * v == F * v
+    @test FCiip * v == FCM * v == F * v == FCMv * v
+    @test FCM2v * v == F * F * v == FCM2iipv * v
     @test @inferred (F * F) * v == @inferred F * (F * v)
     @test @inferred (F * A) * v == @inferred F * (A * v)
     @test LinearMaps._compositemul!(zero(F * A * v), F * A, v, zero(A*v)) ≈ (F * A) * v
@@ -123,6 +128,7 @@ using LinearMaps: LinearMapVector, LinearMapTuple
     w1 = im.*ones(ComplexF64, prod(sizes[1]))
     for i = N:-1:1
         v2 = prod(Lf[i:N]) * ones(prod(sizes[1]))
+        i < N && (y2 = LinearMaps._compositemul!(zero(v2), prod(Lf[i:N]), ones(prod(sizes[1]))))
         u2 = transpose(LinearMap(prod(Lt[N:-1:i]))) * ones(prod(sizes[1]))
         w2 = adjoint(LinearMap(prod(Lc[N:-1:i]))) * ones(prod(sizes[1]))
 
@@ -131,6 +137,7 @@ using LinearMaps: LinearMapVector, LinearMapTuple
         w1 = adjoint(Lc[i]) * w1
 
         @test v1 == v2
+        i < N && @test v2 == y2
         @test u1 == u2
         @test w1 == w2
     end