improved TTFX in multiplication (#164)

Author: jishnub

src/LinearAlgebra/standardchop.jl

@@ -57,8 +57,8 @@ function standardchoplength(coeffs, tol)
     # Step 1: Convert COEFFS to a new monotonically nonincreasing
     #         vector ENVELOPE normalized to begin with the value 1.
 
-    b = abs.(coeffs)
-    m = b[end]*fill(1.0,n)
+    b = map(abs, coeffs)
+    m = fill(b[end]*1.0, n)
     for j = n-1:-1:1
         m[j] = max(b[j], m[j+1]);
     end
@@ -94,13 +94,16 @@ function standardchoplength(coeffs, tol)
         return plateauPoint
     end
 
-    j3 = sum(envelope .≥ tol^(7/6))
+    j3 = count(≥(tol^(7/6)), envelope)
     if j3 < j2
         j2 = j3 + 1
         envelope[j2] = tol^(7/6)
     end
-    cc = log10.(envelope[1:j2])
-    cc .+= range(0, stop=(-1/3)*log10(tol), length=j2)
+    cc = map(log10, @view envelope[1:j2])
+    rng = range(0, stop=(-1/3)*log10(tol), length=j2)
+    for (i,ri) in enumerate(rng)
+        cc[i] += ri
+    end
     d = argmin(cc)
     return max(d - 1, 1)
 end

src/Operators/general/algebra.jl

@@ -63,6 +63,7 @@ domain(P::PlusOperator) = commondomain(P.ops)
 
 _promote_eltypeof(As...) = _promote_eltypeof(As)
 _promote_eltypeof(As::Union{Vector, Tuple}) = mapreduce(eltype, promote_type, As)
+_promote_eltypeof(As::Vector{Operator{T}}) where {T} = T
 
 _extractops(A, ::Any) = [A]
 _extractops(A::PlusOperator, ::typeof(+)) = A.ops
@@ -265,38 +266,26 @@ end
 
 
 
-function promotetimes(opsin::Vector{B},dsp) where B<:Operator
+function promotetimes(opsin::Vector{<:Operator}, dsp = domainspace(last(opsin)))
+
+    @assert length(opsin) > 1 "need at least 2 operators"
     ops=Vector{Operator{_promote_eltypeof(opsin)}}(undef,0)
+    sizehint!(ops, length(opsin))
 
     for k=length(opsin):-1:1
         if !isa(opsin[k],Conversion)
             op=promotedomainspace(opsin[k],dsp)
-            if op==()
-                # do nothing
-            elseif isa(op,TimesOperator)
-                for j=length(op.ops):-1:1
-                    push!(ops,op.ops[j])
-                end
-                dsp=rangespace(op)
+            dsp=rangespace(op)
+            if isa(op,TimesOperator)
+                append!(ops, view(op.ops, reverse(axes(op.ops,1))))
             else
                 push!(ops,op)
-                dsp=rangespace(op)
             end
         end
     end
-    if isempty(ops)
-        ConstantOperator(1.0,dsp)
-    elseif length(ops)==1
-        first(ops)
-    else
-        TimesOperator(reverse!(ops))  # change order in TImesOperator if this is slow
-    end
+    TimesOperator(reverse!(ops))
 end
 
-promotetimes(opsin::Vector{B}) where {B<:Operator}=promotetimes(opsin,domainspace(last(opsin)))
-
-
-
 domainspace(P::TimesOperator)=domainspace(last(P.ops))
 rangespace(P::TimesOperator)=rangespace(first(P.ops))
 
@@ -490,7 +479,8 @@ function *(A::Operator,B::Operator)
     elseif isconstop(B)
         promotedomainspace(strictconvert(Number,B)*A,domainspace(B))
     else
-        promotetimes(Operator{_promote_eltypeof(A, B)}[_extractops(A, *); _extractops(B, *)])
+        promotetimes([_extractops(A, *); _extractops(B, *)],
+            domainspace(B))
     end
 end
 
@@ -619,7 +609,7 @@ function promotedomainspace(P::TimesOperator,sp::Space,cursp::Space)
     elseif length(P.ops)==2
         P.ops[1]*promotedomainspace(P.ops[end],sp)
     else
-        promotetimes([P.ops[1:end-1];promotedomainspace(P.ops[end],sp)])
+        promotetimes([P.ops[1:end-1];promotedomainspace(P.ops[end],sp)], sp)
     end
 end