Avoid some allocations in mutating methods (#347)

* WIP: proof of concept: avoid allocations in inplace methods

* Move zero! methods from arithmetic.jl to functions.jl

* Substitute zero! methods where appropriate

* Substitute 0:a.order -> eachindex(a)

* Add one-order zero! methods

* Try allocation-free zero! for non-mixtures

* WIP: debugging zero stuff

* WIP: workaround TaylorIntegration bug

* Add tests

* Fix typo

* Small fix in div!

* Allocate less in div!(::Taylor1,::NumberNotSeries,::Taylor1,Int)

* Add missing div! method

* Add /(::NumberNotSeries,::Taylor1{TaylorN{<:NumberNotSeries}}

* Update comments

* Address review by @lbenet

* WIP: try less allocations in div!

* Fixes in division methods

* Cleanup and add test

* Cleanup

* More cleanup

* Remove comments

* Update Project.toml: bump minor version and use hyphen instead of commas in IntervalArithmetic compat entry

Project.toml

@@ -1,27 +1,30 @@
 name = "TaylorSeries"
 uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
-version = "0.16.0"
+version = "0.17.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[weakdeps]
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+
+[extensions]
+TaylorSeriesIAExt = "IntervalArithmetic"
+
 [compat]
 Aqua = "0.8"
-IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
+IntervalArithmetic = "0.15 - 0.20"
 LinearAlgebra = "<0.0.1, 1"
 Markdown = "<0.0.1, 1"
 Requires = "0.5.2, 1"
 SparseArrays = "<0.0.1, 1"
 Test = "<0.0.1, 1"
 julia = "1"
 
-[extensions]
-TaylorSeriesIAExt = "IntervalArithmetic"
-
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
@@ -31,6 +34,3 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["IntervalArithmetic", "LinearAlgebra", "SparseArrays", "Test", "Aqua"]
-
-[weakdeps]
-IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"

src/arithmetic.jl

@@ -482,7 +482,7 @@ end
 function mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}},
         ordT::Int) where {T<:NumberNotSeries}
     # Sanity
-    zero!(res, a, ordT)
+    zero!(res, ordT)
     for k in 0:ordT
         @inbounds for ordQ in eachindex(a[ordT])
             mul!(res[ordT], a[k], b[ordT-k], ordQ)
@@ -491,6 +491,27 @@ function mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}}, b::Taylor1{Taylo
     return nothing
 end
 
+@inline function mul!(res::Taylor1{TaylorN{T}}, a::NumberNotSeries,
+    b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
+for l in eachindex(b[k])
+    for m in eachindex(b[k][l])
+        res[k][l][m] = a*b[k][l][m]
+    end
+end
+return nothing
+end
+
+mul!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+b::NumberNotSeries, k::Int) where {T<:NumberNotSeries} = mul!(res, b, a, k)
+
+# in-place division (assumes equal order among TaylorNs)
+function mul!(c::TaylorN, a::TaylorN, b::TaylorN)
+    for k in eachindex(c)
+        mul!(c, a, b, k)
+    end
+end
+
+
 
 @doc doc"""
     mul!(c, a, b, k::Int) --> nothing
@@ -665,6 +686,27 @@ function /(a::Taylor1{TaylorN{T}}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSe
     return res
 end
 
+function /(a::S, b::Taylor1{TaylorN{T}}) where {S<:NumberNotSeries, T<:NumberNotSeries}
+    R = promote_type(TaylorN{S}, TaylorN{T})
+    res = convert(Taylor1{R}, zero(b))
+    iszero(a) && !iszero(b) && return res
+
+    for ordT in eachindex(res)
+        div!(res, a, b, ordT)
+    end
+    return res
+end
+
+function /(a::TaylorN{T}, b::Taylor1{TaylorN{T}}) where {T<:NumberNotSeries}
+    res = zero(b)
+    iszero(a) && !iszero(b) && return res
+
+    aa = Taylor1(a, b.order)
+    for ordT in eachindex(res)
+        div!(res, aa, b, ordT)
+    end
+    return res
+end
 
 @inline function divfactorization(a1::Taylor1, b1::Taylor1)
     # order of first factorized term; a1 and b1 assumed to be of the same order
@@ -731,23 +773,94 @@ end
     return nothing
 end
 
-div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) =
-    div!(v::Taylor1, Taylor1(b, a.order), a, k)
+@inline function div!(c::Taylor1{T}, a::NumberNotSeries,
+        b::Taylor1{T}, k::Int) where {T<:Number}
+    iszero(a) && !iszero(b) && zero!(c, k)
+    # order and coefficient of first factorized term
+    # In this case, since a[k]=0 for k>0, we can simplify to:
+    # ordfact, cdivfact = 0, a/b[0]
+    if k == 0
+        @inbounds c[0] = a/b[0]
+        return nothing
+    end
+
+    @inbounds c[k] = c[0] * b[k]
+    @inbounds for i = 1:k-1
+        c[k] += c[i] * b[k-i]
+    end
+    @inbounds c[k] = -c[k]/b[0]
+    return nothing
+end
+
+# TODO: get rid of remaining allocations here.
+#       This can either be achieved via pre-allocating auxiliary variables
+#       with can be passed here as arguments, or adding allocation-free in-place
+#       methods for operations such as `a += a*b` and `a = -a/b`, where a and b are
+#       TaylorN variables. See #347 for further discussion.
+@inline function div!(c::Taylor1{TaylorN{T}}, a::NumberNotSeries,
+        b::Taylor1{TaylorN{T}}, k::Int) where {T<:NumberNotSeries}
+    iszero(a) && !iszero(b) && zero!(c, k)
+    # order and coefficient of first factorized term
+    # In this case, since a[k]=0 for k>0, we can simplify to:
+    # ordfact, cdivfact = 0, a/b[0]
+    if k == 0
+        @inbounds div!(c[0], a, b[0])
+        return nothing
+    end
+
+    @inbounds mul!(c[k], c[0], b[k])
+    @inbounds for i = 1:k-1
+        c[k] += c[i] * b[k-i]
+    end
+    @inbounds c[k] = -c[k]/b[0]
+    return nothing
+end
 
 # NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0)
 @inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)
     if k==0
-        @inbounds c[0] = a[0] / constant_term(b)
+        @inbounds c[0][1] = constant_term(a) / constant_term(b)
         return nothing
     end
 
     @inbounds for i = 0:k-1
         mul!(c[k], c[i], b[k-i])
     end
-    @inbounds c[k] = (a[k] - c[k]) / constant_term(b)
+    @inbounds for i in eachindex(c[k])
+        c[k][i] = (a[k][i] - c[k][i]) / constant_term(b)
+    end
     return nothing
 end
 
+# NOTE: Here `div!` *accumulates* the result of a / b in c[k] (k > 0)
+@inline function div!(c::TaylorN, a::NumberNotSeries, b::TaylorN, k::Int)
+    if k==0
+        @inbounds c[0][1] = a / constant_term(b)
+        return nothing
+    end
+
+    @inbounds for i = 0:k-1
+        mul!(c[k], c[i], b[k-i])
+    end
+    @inbounds for i in eachindex(c[k])
+        c[k][i] = ( -c[k][i] ) / constant_term(b)
+    end
+    return nothing
+end
+
+# in-place division (assumes equal order among TaylorNs)
+function div!(c::TaylorN, a::TaylorN, b::TaylorN)
+    for k in eachindex(c)
+        div!(c, a, b, k)
+    end
+end
+
+function div!(c::TaylorN, a::NumberNotSeries, b::TaylorN)
+    for k in eachindex(c)
+        div!(c, a, b, k)
+    end
+end
+
 # NOTE: Here `div!` *accumulates* the result of a / b in res[k] (k > 0)
 @inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
         b::Taylor1{TaylorN{T}}, ordT::Int) where {T<:NumberNotSeriesN}
@@ -758,7 +871,7 @@ end
         @inbounds res[0] = cdivfact
         return nothing
     end
-    zero!(res, a, ordT)
+    zero!(res, ordT)
     imin = max(0, ordT+ordfact-b.order)
     aux = TaylorN(zero(a[ordT][0][1]), a[ordT].order)
     for k in imin:ordT-1
@@ -783,6 +896,15 @@ end
     return nothing
 end
 
+@inline function div!(res::Taylor1{TaylorN{T}}, a::Taylor1{TaylorN{T}},
+        b::NumberNotSeries, k::Int) where {T<:NumberNotSeries}
+    for l in eachindex(a[k])
+        for m in eachindex(a[k][l])
+            res[k][l][m] = a[k][l][m]/b
+        end
+    end
+    return nothing
+end