Fix Jacobi deriv evaluate for endpoints

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.5.14"
+version = "0.5.15"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/Spaces/Chebyshev/ChebyshevOperators.jl

@@ -17,9 +17,13 @@ Evaluation(S::Chebyshev,x::typeof(leftendpoint),o::Integer) =
     ConcreteEvaluation(S,x,o)
 Evaluation(S::Chebyshev,x::typeof(rightendpoint),o::Integer) =
     ConcreteEvaluation(S,x,o)
+Evaluation(S::Chebyshev,x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
 
-Evaluation(S::Chebyshev,x::Number,o::Integer) =
-    o==0 ? ConcreteEvaluation(S,x,o) : EvaluationWrapper(S,x,o,Evaluation(x)*Derivative(S,o))
+Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::typeof(leftendpoint),o::Integer) =
+    ConcreteEvaluation(S,x,o)
+Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::typeof(rightendpoint),o::Integer) =
+    ConcreteEvaluation(S,x,o)
+Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
 
 function evaluatechebyshev(n::Integer,x::T) where T<:Number
     if n == 1
@@ -28,15 +32,22 @@ function evaluatechebyshev(n::Integer,x::T) where T<:Number
         p = zeros(T,n)
         p[1] = one(T)
         p[2] = x
+        twox = 2x
 
         for j=2:n-1
-            p[j+1] = 2x*p[j] - p[j-1]
+            p[j+1] = muladd(twox, p[j], -p[j-1])
         end
 
         p
     end
 end
 
+# We assume that x is already scaled to the canonical domain. S is unused here
+function forwardrecurrence(::Type{T},S::Chebyshev,r::AbstractUnitRange{<:Integer},x::Number) where {T}
+    @assert !isempty(r) && first(r) >= 0
+    v = evaluatechebyshev(maximum(r)+1, T(x))
+    first(r) == 0 ? v : v[r .+ 1]
+end
 
 function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
     T=eltype(op)
@@ -60,7 +71,8 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint
     end
 end
 
-function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
+function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},k::AbstractUnitRange) where {DD<:IntervalOrSegment,RR}
+    Base.require_one_based_indexing(k)
     T=eltype(op)
     x = op.x
     d = domain(op)
@@ -69,15 +81,13 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)
     n=length(k)
 
     ret = Array{T}(undef, n)
-    k1=1-first(k)
-    @simd for j=k
-        @inbounds ret[j+k1]=(-1)^(p+1)*(-one(T))^j
+    for (ind, j) in enumerate(k)
+        ret[ind]=(-1)^(p+1)*(-one(T))^j
     end
 
-    for m=0:p-1
-        k1=1-first(k)
-        @simd for j=k
-            @inbounds ret[j+k1] *= (j-1)^2-m^2
+    for m in 0:p-1
+        for (ind, j) in enumerate(k)
+            ret[ind] *= (j-1)^2-m^2
         end
         scal!(strictconvert(T,1/(2m+1)), ret)
     end
@@ -86,7 +96,8 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)
 end
 
 
-function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
+function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},k::AbstractUnitRange) where {DD<:IntervalOrSegment,RR}
+    Base.require_one_based_indexing(k)
     T=eltype(op)
     x = op.x
     d = domain(op)
@@ -96,35 +107,16 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint
 
     ret = fill(one(T),n)
 
-    for m=0:p-1
-        k1=1-first(k)
-        @simd for j=k
-            @inbounds ret[j+k1] *= (j-1)^2-m^2
+    for m in 0:p-1
+        for (ind, j) in enumerate(k)
+            ret[ind] *= (j-1)^2-m^2
         end
         scal!(strictconvert(T,1/(2m+1)), ret)
     end
 
     scal!(cst,ret)
 end
 
-function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},M,OT,T},
-                j::Integer) where {DD<:IntervalOrSegment,RR,M<:Real,OT,T}
-    if op.order == 0
-        strictconvert(T,evaluatechebyshev(j,tocanonical(domain(op),op.x))[end])
-    else
-        error("Only zero–second order implemented")
-    end
-end
-
-function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},M,OT,T},
-                                             k::AbstractRange) where {DD<:IntervalOrSegment,RR,M<:Real,OT,T}
-    if op.order == 0
-        Array{T}(evaluatechebyshev(k[end],tocanonical(domain(op),op.x))[k])
-    else
-        error("Only zero–second order implemented")
-    end
-end
-
 @inline function _Dirichlet_Chebyshev(S, order)
     order == 0 && return ConcreteDirichlet(S,ArraySpace([ConstantSpace.(Point.(endpoints(domain(S))))...]),0)
     default_Dirichlet(S,order)

src/Spaces/Jacobi/JacobiOperators.jl

@@ -1,99 +1,3 @@
-## Evaluation
-
-
-function Evaluation(S::Jacobi,x::Union{typeof(leftendpoint),typeof(leftendpoint)},order)
-    if order ≤ 2
-        ConcreteEvaluation(S,x,order)
-    else
-        # assume Derivative is available
-        D = Derivative(S,order)
-        EvaluationWrapper(S,x,order,Evaluation(rangespace(D),x)*D)
-    end
-end
-
-function Evaluation(S::Jacobi,x,order)
-    if order == 0
-        ConcreteEvaluation(S,x,order)
-    else
-        # assume Derivative is available
-        D = Derivative(S,order)
-        EvaluationWrapper(S,x,order,Evaluation(rangespace(D),x)*D)
-    end
-end
-
-# avoid ambiguity
-for OP in (:leftendpoint,:rightendpoint)
-    @eval getindex(op::ConcreteEvaluation{<:Jacobi,typeof($OP)},k::Integer) =
-        op[k:k][1]
-end
-
-getindex(op::ConcreteEvaluation{<:Jacobi},k::Integer) = op[k:k][1]
-
-
-function getindex(op::ConcreteEvaluation{<:Jacobi,typeof(leftendpoint)},kr::AbstractRange)
-    @assert op.order <= 2
-    sp=op.space
-    T=eltype(op)
-    RT=real(T)
-    a=strictconvert(RT,sp.a);b=strictconvert(RT,sp.b)
-
-    if op.order == 0
-        jacobip(T,kr.-1,a,b,-one(T))
-    elseif op.order == 1&&  b==0
-        d=domain(op)
-        @assert isa(d,IntervalOrSegment)
-        T[tocanonicalD(d,leftendpoint(d))/2*(a+k)*(k-1)*(-1)^k for k=kr]
-    elseif op.order == 1
-        d=domain(op)
-        @assert isa(d,IntervalOrSegment)
-        if kr[1]==1 && kr[end] ≥ 2
-            tocanonicalD(d,leftendpoint(d)).*(a .+ b .+ kr).*T[zero(T);jacobip(T,0:kr[end]-2,1+a,1+b,-one(T))]/2
-        elseif kr[1]==1  # kr[end] ≤ 1
-            zeros(T,length(kr))
-        else
-            tocanonicalD(d,leftendpoint(d)) .* (a.+b.+kr).*jacobip(T,kr.-1,1+a,1+b,-one(T))/2
-        end
-    elseif op.order == 2
-        @assert b==0
-        @assert domain(op) == ChebyshevInterval()
-        T[-0.125*(a+k)*(a+k+1)*(k-2)*(k-1)*(-1)^k for k=kr]
-    else
-        error("Not implemented")
-    end
-end
-
-function getindex(op::ConcreteEvaluation{<:Jacobi,typeof(rightendpoint)},kr::AbstractRange)
-    @assert op.order <= 2
-    sp=op.space
-    T=eltype(op)
-    RT=real(T)
-    a=strictconvert(RT,sp.a);b=strictconvert(RT,sp.b)
-
-
-    if op.order == 0
-        jacobip(T,kr.-1,a,b,one(T))
-    elseif op.order == 1
-        d=domain(op)
-        @assert isa(d,IntervalOrSegment)
-        if kr[1]==1 && kr[end] ≥ 2
-            tocanonicalD(d,leftendpoint(d))*((a+b).+kr).*T[zero(T);jacobip(T,0:kr[end]-2,1+a,1+b,one(T))]/2
-        elseif kr[1]==1  # kr[end] ≤ 1
-            zeros(T,length(kr))
-        else
-            tocanonicalD(d,leftendpoint(d))*((a+b).+kr).*jacobip(T,kr.-1,1+a,1+b,one(T))/2
-        end
-    else
-        error("Not implemented")
-    end
-end
-
-
-function getindex(op::ConcreteEvaluation{<:Jacobi,<:Number},kr::AbstractRange)
-    @assert op.order == 0
-    jacobip(eltype(op),kr.-1,op.space.a,op.space.b,tocanonical(domain(op),op.x))
-end
-
-
 ## Derivative
 
 function Derivative(J::Jacobi,k::Integer)
@@ -113,7 +17,19 @@ getindex(T::ConcreteDerivative{J},k::Integer,j::Integer) where {J<:Jacobi} =
     j==k+1 ? eltype(T)((k+1+T.space.a+T.space.b)/complexlength(domain(T))) : zero(eltype(T))
 
 
+# Evaluation
+
+Evaluation(S::Jacobi,x::typeof(leftendpoint),o::Integer) =
+    ConcreteEvaluation(S,x,o)
+Evaluation(S::Jacobi,x::typeof(rightendpoint),o::Integer) =
+    ConcreteEvaluation(S,x,o)
+Evaluation(S::Jacobi,x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
 
+Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::typeof(leftendpoint),o::Integer) =
+    ConcreteEvaluation(S,x,o)
+Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::typeof(rightendpoint),o::Integer) =
+    ConcreteEvaluation(S,x,o)
+Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
 
 ## Integral
 

src/Spaces/PolynomialSpace.jl

@@ -301,39 +301,74 @@ function forwardrecurrence(::Type{T},S::Space,r::AbstractRange,x::Number) where
     return v[r.+1]
 end
 
-
-function Evaluation(S::PolynomialSpace,x,order)
-    if order == 0
-        ConcreteEvaluation(S,x,order)
-    else
-        # assume Derivative is available
-        D = Derivative(S,order)
-        EvaluationWrapper(S,x,order,Evaluation(rangespace(D),x)*D)
-    end
+getindex(op::ConcreteEvaluation{<:PolynomialSpace}, k::Integer) = op[k:k][1]
+# avoid ambiguity
+for OP in (:leftendpoint, :rightendpoint)
+    @eval getindex(op::ConcreteEvaluation{<:PolynomialSpace,typeof($OP)},k::Integer) =
+        op[k:k][1]
 end
 
-
 function getindex(op::ConcreteEvaluation{J,typeof(leftendpoint)},kr::AbstractRange) where J<:PolynomialSpace
-    sp=op.space
-    T=eltype(op)
-    @assert op.order == 0
-    forwardrecurrence(T,sp,kr.-1,-one(T))
+    _getindex_evaluation(op, leftendpoint(domain(op)), kr)
 end
 
 function getindex(op::ConcreteEvaluation{J,typeof(rightendpoint)},kr::AbstractRange) where J<:PolynomialSpace
-    sp=op.space
-    T=eltype(op)
-    @assert op.order == 0
-    forwardrecurrence(T,sp,kr.-1,one(T))
+    _getindex_evaluation(op, rightendpoint(domain(op)), kr)
 end
 
+function getindex(op::ConcreteEvaluation{J,TT}, kr::AbstractRange) where {J<:PolynomialSpace,TT<:Number}
+    _getindex_evaluation(op, op.x, kr)
+end
+
+function _getindex_evaluation(op::ConcreteEvaluation{<:PolynomialSpace}, x, kr::AbstractRange)
+    _getindex_evaluation(eltype(op), op.space, op.order, x, kr)
+end
 
-function getindex(op::ConcreteEvaluation{J,TT},kr::AbstractRange) where {J<:PolynomialSpace,TT<:Number}
-    sp=op.space
-    T=eltype(op)
-    x=op.x
-    @assert op.order == 0
-    forwardrecurrence(T,sp,kr .- 1,tocanonical(sp,x))
+function _getindex_evaluation(::Type{T}, sp, order, x, kr::AbstractRange) where {T}
+    Base.require_one_based_indexing(kr)
+    isempty(kr) && return zeros(T, 0)
+    if order == 0
+        forwardrecurrence(T,sp,kr .- 1,tocanonical(sp,x))
+    else
+        z = Zeros{T}(length(range(minimum(kr), order, step=step(kr))))
+        kr_red = kr .- (order + 1)
+        polydegrees = reverse(range(maximum(kr_red), max(0, minimum(kr_red)), step=-step(kr)))
+        D = Derivative(sp, order)
+        if !isempty(polydegrees)
+            P = forwardrecurrence(T, rangespace(D), polydegrees, tocanonical(sp,x))
+            # in case the derivative has only one non-zero band, the matrix-vector
+            # multiplication simplifies considerably.
+            # This branch is particularly useful for ConcreteDerivatives where
+            # indexing is fast, as the non-zero band may be computed without
+            # evaluating the matrix
+            if bandwidth(D, 1) == -bandwidth(D, 2)
+                d = nonzeroband(T, D, polydegrees, order)
+                Pd = P .* d
+                if !isempty(z)
+                    Pd = T[z; Pd]
+                end
+            else
+                # in general, this is a banded matrix-vector product
+                Dtv = view(transpose(D), kr, 1:length(P))
+                Pd = strictconvert(Vector{T},  mul_coefficients(Dtv, P))
+            end
+        else
+            Pd = T[z;]
+        end
+        Pd
+    end
+end
+
+function nonzeroband(::Type{T}, D::ConcreteDerivative, polydegrees, order) where {T}
+    bw = bandwidth(D, 2)
+    T[D[k+1, k+1+bw] for k in polydegrees]
+end
+function nonzeroband(::Type{T}, D, polydegrees, order) where {T}
+    bw = bandwidth(D, 2)
+    rows = 1:(maximum(polydegrees) + order + 1)
+    B = D[rows, rows .+ bw]
+    Bv = @view B[diagind(B)]
+    Bv[polydegrees .+ 1]
 end