lincomb optimization, added some TODOs in light of mulstyle

Author: dkarrasch

src/linearcombination.jl

@@ -66,23 +66,35 @@ LinearAlgebra.adjoint(A::LinearCombination)   = LinearCombination{eltype(A)}(map
 # map types that have an allocation-free 5-arg mul! implementation
 const FreeMap = Union{MatrixMap,UniformScalingMap}
 
+# with the mulstyle feature, this version distinction becomes obsolete
+# the function for older versions would automatically call the ThreeArg method
+# of _lincombmul!
 if VERSION < v"1.3.0-alpha.115"
 
 Base.@propagate_inbounds function LinearAlgebra.mul!(y::AbstractVector, A::LinearCombination, x::AbstractVector,
                     α::Number=true, β::Number=false)
     @boundscheck check_dim_mul(y, A, x)
-    mul!(y, A.maps[1], x, α, β)
-    @inbounds begin
-        l = length(A.maps)
-        if l>1
-            z = similar(y)
-            for n in 2:l
-                mul!(z, A.maps[n], x, α, false)
-                y .+= z
+    if iszero(α)
+        iszero(β) && (fill!(y, zero(eltype(y))); return y)
+        isone(β) && return y
+        # β != 0, 1
+        rmul!(y, β)
+        return y
+    else
+        # TODO: reconsider the below code in light of the mulstyle trait
+        mul!(y, A.maps[1], x, α, β)
+        @inbounds begin
+            l = length(A.maps)
+            if l>1
+                z = similar(y)
+                for n in 2:l
+                    mul!(z, A.maps[n], x, α, false)
+                    y .+= z
+                end
             end
         end
+        return y
     end
-    return y
 end
 
 else # 5-arg mul! is available for matrices
@@ -98,30 +110,58 @@ end # VERSION
 ############
 # multiplication helper functions
 ############
-
+# to be replaced by _lincombmul!(y, A::LinearCombination{<:Any,FiveArg}, x, α::Number, β::Number)
 @inline function _lincombmul!(y, A::LinearCombination{<:Any,<:Tuple{Vararg{FreeMap}}}, x,
                                 α::Number, β::Number)
-    mul!(y, A.maps[1], x, α, β)
-    @inbounds for n in 2:length(A.maps)
-        mul!(y, A.maps[n], x, α, true)
+    if iszero(α)
+        iszero(β) && (fill!(y, zero(eltype(y))); return y)
+        isone(β) && return y
+        # β != 0, 1
+        rmul!(y, β)
+        return y
+    else
+        mul!(y, A.maps[1], x, α, β)
+        @inbounds for n in 2:length(A.maps)
+            mul!(y, A.maps[n], x, α, true)
+        end
+        return y
     end
-    return y
 end
 
+# to be replaced by _lincombmul!(y, A::LinearCombination{<:Any,ThreeArg}, x, α::Number, β::Number)
 @inline function _lincombmul!(y, A::LinearCombination, x, α::Number, β::Number)
-    mul!(y, first(A.maps), x, α, β)
-    l = length(A.maps)
-    if l>1
+    if iszero(α)
+        iszero(β) && (fill!(y, zero(eltype(y))); return y)
+        isone(β) && return y
+        # β != 0, 1
+        rmul!(y, β)
+        return y
+    else
         z = similar(y)
-        for n in 2:l
-            @inbounds An = A.maps[n]
-            if An isa FreeMap
-                mul!(y, An, x, α, true)
+        A1 = A.maps[1]
+        if A1 isa FreeMap # to be replaced by mulstyle(A1) == FiveArg
+            mul!(y, A1, x, α, β)
+        else
+            if iszero(β)
+                mul!(y, A1, x)
+                !isone(α) && rmul!(y, α)
             else
-                mul!(z, An, x, α, false)
-                y .+= z
+                !isone(β) && rmul!(y, β)
+                mul!(z, A1, x)
+                y .+= z .* α
+        end
+
+        l = length(A.maps)
+        if l>1
+            @inbounds for An in Base.tail(A.maps)
+                if An isa FreeMap # to be replaced by mulstyle(An) == FiveArg
+                    mul!(y, An, x, α, true)
+                else
+                    mul!(z, An, x)
+                    y .+= z .* α
+                end
             end
         end
+        return y
     end
-    return y
 end