reduce allocations in Chebyshev deriv (#280)

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.39"
+version = "0.6.40"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"

src/Spaces/Chebyshev/Chebyshev.jl

@@ -204,17 +204,18 @@ end
 
 
 # diff T -> U, then convert U -> T
-function integrate(f::Fun{Chebyshev{D,R}}) where {D<:IntervalOrSegment,R}
+function integrate(f::Fun{<:Chebyshev{<:IntervalOrSegment}})
     cfs = coefficients(f)
     z = fromcanonicalD(f,0)
     v = z .* float.(cfs)
     ultraint!(ultraconversion!(v))
     Fun(f.space, v)
 end
-function differentiate(f::Fun{Chebyshev{D,R}}) where {D<:IntervalOrSegment,R}
+function differentiate(f::Fun{<:Chebyshev{<:IntervalOrSegment}})
     cfs = coefficients(f)
     z = fromcanonicalD(f,0)
-    v = ultraiconversion(ultradiff(cfs)) ./ z
+    w = float.(cfs) ./ z
+    v = ultraiconversion!(ultradiff!(w))
     Fun(f.space, v)
 end
 

src/ultraspherical.jl

@@ -3,19 +3,22 @@ export ultraconversion!,ultraint!
 ## Start of support for UFun
 
 # diff from T -> U
-function ultradiff(v::AbstractVector{T}) where T<:Number
+function ultradiff(v::AbstractVector{<:Number})
+    ultradiff!(copy(v))
+end
+
+function ultradiff!(v::AbstractVector{<:Number})
     Base.require_one_based_indexing(v)
     #polynomial is p(x) = sum ( v[i] * x^(i-1) )
-    if length(v)≤1
-        w = zeros(T,1)
-    else
-        w = Array{T}(undef, length(v)-1)
-        for k in eachindex(w)
-            @inbounds w[k] = k*v[k+1]
-        end
+    if length(v) <= 1
+        fill!(v, zero(eltype(v)))
+        return v
     end
-
-    w
+    for k in eachindex(v)[1:end-1]
+        v[k] = k*v[k+1]
+    end
+    resize!(v, length(v)-1)
+    return v
 end
 
 #int from U ->T
@@ -47,31 +50,30 @@ function ultraint!(v::AbstractVector{T}) where T<:Number
 end
 
 # Convert from U -> T
-function ultraiconversion(v::AbstractVector{T}) where T<:Number
+function ultraiconversion(v::AbstractVector{<:Number})
+    ultraiconversion!(copy(v))
+end
+
+function ultraiconversion!(v::AbstractVector{<:Number})
     Base.require_one_based_indexing(v)
     n = length(v)
-    w = Array{T}(undef, n)
 
-    if n == 1
-        w[1] = v[1]
-    elseif n == 2
-        w[1] = v[1]
-        w[2] = 2v[2]
+    if n == 2
+        @inbounds v[2] *= 2
     elseif n ≥ 3
-        @inbounds w[end] = 2v[end]
-        @inbounds w[end-1] = 2v[end-1]
+        @inbounds v[end] *= 2
+        @inbounds v[end-1] *= 2
 
         for k = n-2:-1:2
-            @inbounds w[k] = 2*(v[k] + .5w[k+2])
+            @inbounds v[k] = 2v[k] + v[k+2]
         end
 
-        @inbounds w[1] = v[1] + .5w[3]
+        @inbounds v[1] += 0.5 * v[3]
     end
 
-    w
+    return v
 end
 
-
 # Convert T -> U
 function ultraconversion(v::AbstractVector{<:Number})
     ultraconversion!(float.(v))

test/ChebyshevTest.jl

@@ -117,6 +117,12 @@ include("testutils.jl")
 
         @test @inferred(sum(ef))  ≈ 2.3504023872876028
         @test @inferred(norm(ef))  ≈ 1.90443178083307
+
+        @testset "edge cases" begin
+            @test Fun(Chebyshev(), [1.0])' ≈ 0
+            @test Fun(Chebyshev(), Float64[])' ≈ 0
+            @test Fun(x->x^3, Chebyshev(2..5))' ≈ Fun(x->3x^2, Chebyshev(2..5))
+        end
     end
 
     @testset "other domains" begin