Skip to content

Commit 0cf2115

Browse files
jishnubjishnub
authored
Fix Jacobi deriv evaluate for endpoints (#142)
* Fix Jacobi deriv evaluate for endpoints * simplify branches * fix range on v1.6 * concrete evaluate for Polynomialspace * specialize _getindex_evaluation * Fixes for Chebyshev Evaluation * fix tocanonical in forwardrecursion * Approximate check in indexing * indexing for banded derivatives
1 parent 968475f commit 0cf2115

File tree

9 files changed

+266
-161
lines changed

9 files changed

+266
-161
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "ApproxFunOrthogonalPolynomials"
22
uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
3-
version = "0.5.14"
3+
version = "0.5.15"
44

55
[deps]
66
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

‎src/Spaces/Chebyshev/ChebyshevOperators.jl‎

Lines changed: 26 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,13 @@ Evaluation(S::Chebyshev,x::typeof(leftendpoint),o::Integer) =
1717
ConcreteEvaluation(S,x,o)
1818
Evaluation(S::Chebyshev,x::typeof(rightendpoint),o::Integer) =
1919
ConcreteEvaluation(S,x,o)
20+
Evaluation(S::Chebyshev,x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
2021

21-
Evaluation(S::Chebyshev,x::Number,o::Integer) =
22-
o==0 ? ConcreteEvaluation(S,x,o) : EvaluationWrapper(S,x,o,Evaluation(x)*Derivative(S,o))
22+
Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::typeof(leftendpoint),o::Integer) =
23+
ConcreteEvaluation(S,x,o)
24+
Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::typeof(rightendpoint),o::Integer) =
25+
ConcreteEvaluation(S,x,o)
26+
Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
2327

2428
function evaluatechebyshev(n::Integer,x::T) where T<:Number
2529
if n == 1
@@ -28,15 +32,22 @@ function evaluatechebyshev(n::Integer,x::T) where T<:Number
2832
p = zeros(T,n)
2933
p[1] = one(T)
3034
p[2] = x
35+
twox = 2x
3136

3237
for j=2:n-1
33-
p[j+1] = 2x*p[j] - p[j-1]
38+
p[j+1] = muladd(twox, p[j], -p[j-1])
3439
end
3540

3641
p
3742
end
3843
end
3944

45+
# We assume that x is already scaled to the canonical domain. S is unused here
46+
function forwardrecurrence(::Type{T},S::Chebyshev,r::AbstractUnitRange{<:Integer},x::Number) where {T}
47+
@assert !isempty(r) && first(r) >= 0
48+
v = evaluatechebyshev(maximum(r)+1, T(x))
49+
first(r) == 0 ? v : v[r .+ 1]
50+
end
4051

4152
function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
4253
T=eltype(op)
@@ -60,7 +71,8 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint
6071
end
6172
end
6273

63-
function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
74+
function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},k::AbstractUnitRange) where {DD<:IntervalOrSegment,RR}
75+
Base.require_one_based_indexing(k)
6476
T=eltype(op)
6577
x = op.x
6678
d = domain(op)
@@ -69,15 +81,13 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)
6981
n=length(k)
7082

7183
ret = Array{T}(undef, n)
72-
k1=1-first(k)
73-
@simd for j=k
74-
@inbounds ret[j+k1]=(-1)^(p+1)*(-one(T))^j
84+
for (ind, j) in enumerate(k)
85+
ret[ind]=(-1)^(p+1)*(-one(T))^j
7586
end
7687

77-
for m=0:p-1
78-
k1=1-first(k)
79-
@simd for j=k
80-
@inbounds ret[j+k1] *= (j-1)^2-m^2
88+
for m in 0:p-1
89+
for (ind, j) in enumerate(k)
90+
ret[ind] *= (j-1)^2-m^2
8191
end
8292
scal!(strictconvert(T,1/(2m+1)), ret)
8393
end
@@ -86,7 +96,8 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)
8696
end
8797

8898

89-
function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
99+
function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},k::AbstractUnitRange) where {DD<:IntervalOrSegment,RR}
100+
Base.require_one_based_indexing(k)
90101
T=eltype(op)
91102
x = op.x
92103
d = domain(op)
@@ -96,35 +107,16 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint
96107

97108
ret = fill(one(T),n)
98109

99-
for m=0:p-1
100-
k1=1-first(k)
101-
@simd for j=k
102-
@inbounds ret[j+k1] *= (j-1)^2-m^2
110+
for m in 0:p-1
111+
for (ind, j) in enumerate(k)
112+
ret[ind] *= (j-1)^2-m^2
103113
end
104114
scal!(strictconvert(T,1/(2m+1)), ret)
105115
end
106116

107117
scal!(cst,ret)
108118
end
109119

110-
function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},M,OT,T},
111-
j::Integer) where {DD<:IntervalOrSegment,RR,M<:Real,OT,T}
112-
if op.order == 0
113-
strictconvert(T,evaluatechebyshev(j,tocanonical(domain(op),op.x))[end])
114-
else
115-
error("Only zero–second order implemented")
116-
end
117-
end
118-
119-
function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},M,OT,T},
120-
k::AbstractRange) where {DD<:IntervalOrSegment,RR,M<:Real,OT,T}
121-
if op.order == 0
122-
Array{T}(evaluatechebyshev(k[end],tocanonical(domain(op),op.x))[k])
123-
else
124-
error("Only zero–second order implemented")
125-
end
126-
end
127-
128120
@inline function _Dirichlet_Chebyshev(S, order)
129121
order == 0 && return ConcreteDirichlet(S,ArraySpace([ConstantSpace.(Point.(endpoints(domain(S))))...]),0)
130122
default_Dirichlet(S,order)

‎src/Spaces/Jacobi/JacobiOperators.jl‎

Lines changed: 12 additions & 96 deletions
Original file line numberDiff line numberDiff line change
@@ -1,99 +1,3 @@
1-
## Evaluation
2-
3-
4-
function Evaluation(S::Jacobi,x::Union{typeof(leftendpoint),typeof(leftendpoint)},order)
5-
if order ≤ 2
6-
ConcreteEvaluation(S,x,order)
7-
else
8-
# assume Derivative is available
9-
D = Derivative(S,order)
10-
EvaluationWrapper(S,x,order,Evaluation(rangespace(D),x)*D)
11-
end
12-
end
13-
14-
function Evaluation(S::Jacobi,x,order)
15-
if order == 0
16-
ConcreteEvaluation(S,x,order)
17-
else
18-
# assume Derivative is available
19-
D = Derivative(S,order)
20-
EvaluationWrapper(S,x,order,Evaluation(rangespace(D),x)*D)
21-
end
22-
end
23-
24-
# avoid ambiguity
25-
for OP in (:leftendpoint,:rightendpoint)
26-
@eval getindex(op::ConcreteEvaluation{<:Jacobi,typeof($OP)},k::Integer) =
27-
op[k:k][1]
28-
end
29-
30-
getindex(op::ConcreteEvaluation{<:Jacobi},k::Integer) = op[k:k][1]
31-
32-
33-
function getindex(op::ConcreteEvaluation{<:Jacobi,typeof(leftendpoint)},kr::AbstractRange)
34-
@assert op.order <= 2
35-
sp=op.space
36-
T=eltype(op)
37-
RT=real(T)
38-
a=strictconvert(RT,sp.a);b=strictconvert(RT,sp.b)
39-
40-
if op.order == 0
41-
jacobip(T,kr.-1,a,b,-one(T))
42-
elseif op.order == 1&& b==0
43-
d=domain(op)
44-
@assert isa(d,IntervalOrSegment)
45-
T[tocanonicalD(d,leftendpoint(d))/2*(a+k)*(k-1)*(-1)^k for k=kr]
46-
elseif op.order == 1
47-
d=domain(op)
48-
@assert isa(d,IntervalOrSegment)
49-
if kr[1]==1 && kr[end] ≥ 2
50-
tocanonicalD(d,leftendpoint(d)).*(a .+ b .+ kr).*T[zero(T);jacobip(T,0:kr[end]-2,1+a,1+b,-one(T))]/2
51-
elseif kr[1]==1 # kr[end] ≤ 1
52-
zeros(T,length(kr))
53-
else
54-
tocanonicalD(d,leftendpoint(d)) .* (a.+b.+kr).*jacobip(T,kr.-1,1+a,1+b,-one(T))/2
55-
end
56-
elseif op.order == 2
57-
@assert b==0
58-
@assert domain(op) == ChebyshevInterval()
59-
T[-0.125*(a+k)*(a+k+1)*(k-2)*(k-1)*(-1)^k for k=kr]
60-
else
61-
error("Not implemented")
62-
end
63-
end
64-
65-
function getindex(op::ConcreteEvaluation{<:Jacobi,typeof(rightendpoint)},kr::AbstractRange)
66-
@assert op.order <= 2
67-
sp=op.space
68-
T=eltype(op)
69-
RT=real(T)
70-
a=strictconvert(RT,sp.a);b=strictconvert(RT,sp.b)
71-
72-
73-
if op.order == 0
74-
jacobip(T,kr.-1,a,b,one(T))
75-
elseif op.order == 1
76-
d=domain(op)
77-
@assert isa(d,IntervalOrSegment)
78-
if kr[1]==1 && kr[end] ≥ 2
79-
tocanonicalD(d,leftendpoint(d))*((a+b).+kr).*T[zero(T);jacobip(T,0:kr[end]-2,1+a,1+b,one(T))]/2
80-
elseif kr[1]==1 # kr[end] ≤ 1
81-
zeros(T,length(kr))
82-
else
83-
tocanonicalD(d,leftendpoint(d))*((a+b).+kr).*jacobip(T,kr.-1,1+a,1+b,one(T))/2
84-
end
85-
else
86-
error("Not implemented")
87-
end
88-
end
89-
90-
91-
function getindex(op::ConcreteEvaluation{<:Jacobi,<:Number},kr::AbstractRange)
92-
@assert op.order == 0
93-
jacobip(eltype(op),kr.-1,op.space.a,op.space.b,tocanonical(domain(op),op.x))
94-
end
95-
96-
971
## Derivative
982

993
function Derivative(J::Jacobi,k::Integer)
@@ -113,7 +17,19 @@ getindex(T::ConcreteDerivative{J},k::Integer,j::Integer) where {J<:Jacobi} =
11317
j==k+1 ? eltype(T)((k+1+T.space.a+T.space.b)/complexlength(domain(T))) : zero(eltype(T))
11418

11519

20+
# Evaluation
21+
22+
Evaluation(S::Jacobi,x::typeof(leftendpoint),o::Integer) =
23+
ConcreteEvaluation(S,x,o)
24+
Evaluation(S::Jacobi,x::typeof(rightendpoint),o::Integer) =
25+
ConcreteEvaluation(S,x,o)
26+
Evaluation(S::Jacobi,x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
11627

28+
Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::typeof(leftendpoint),o::Integer) =
29+
ConcreteEvaluation(S,x,o)
30+
Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::typeof(rightendpoint),o::Integer) =
31+
ConcreteEvaluation(S,x,o)
32+
Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
11733

11834
## Integral
11935

‎src/Spaces/PolynomialSpace.jl‎

Lines changed: 59 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -301,39 +301,74 @@ function forwardrecurrence(::Type{T},S::Space,r::AbstractRange,x::Number) where
301301
return v[r.+1]
302302
end
303303

304-
305-
function Evaluation(S::PolynomialSpace,x,order)
306-
if order == 0
307-
ConcreteEvaluation(S,x,order)
308-
else
309-
# assume Derivative is available
310-
D = Derivative(S,order)
311-
EvaluationWrapper(S,x,order,Evaluation(rangespace(D),x)*D)
312-
end
304+
getindex(op::ConcreteEvaluation{<:PolynomialSpace}, k::Integer) = op[k:k][1]
305+
# avoid ambiguity
306+
for OP in (:leftendpoint, :rightendpoint)
307+
@eval getindex(op::ConcreteEvaluation{<:PolynomialSpace,typeof($OP)},k::Integer) =
308+
op[k:k][1]
313309
end
314310

315-
316311
function getindex(op::ConcreteEvaluation{J,typeof(leftendpoint)},kr::AbstractRange) where J<:PolynomialSpace
317-
sp=op.space
318-
T=eltype(op)
319-
@assert op.order == 0
320-
forwardrecurrence(T,sp,kr.-1,-one(T))
312+
_getindex_evaluation(op, leftendpoint(domain(op)), kr)
321313
end
322314

323315
function getindex(op::ConcreteEvaluation{J,typeof(rightendpoint)},kr::AbstractRange) where J<:PolynomialSpace
324-
sp=op.space
325-
T=eltype(op)
326-
@assert op.order == 0
327-
forwardrecurrence(T,sp,kr.-1,one(T))
316+
_getindex_evaluation(op, rightendpoint(domain(op)), kr)
328317
end
329318

319+
function getindex(op::ConcreteEvaluation{J,TT}, kr::AbstractRange) where {J<:PolynomialSpace,TT<:Number}
320+
_getindex_evaluation(op, op.x, kr)
321+
end
322+
323+
function _getindex_evaluation(op::ConcreteEvaluation{<:PolynomialSpace}, x, kr::AbstractRange)
324+
_getindex_evaluation(eltype(op), op.space, op.order, x, kr)
325+
end
330326

331-
function getindex(op::ConcreteEvaluation{J,TT},kr::AbstractRange) where {J<:PolynomialSpace,TT<:Number}
332-
sp=op.space
333-
T=eltype(op)
334-
x=op.x
335-
@assert op.order == 0
336-
forwardrecurrence(T,sp,kr .- 1,tocanonical(sp,x))
327+
function _getindex_evaluation(::Type{T}, sp, order, x, kr::AbstractRange) where {T}
328+
Base.require_one_based_indexing(kr)
329+
isempty(kr) && return zeros(T, 0)
330+
if order == 0
331+
forwardrecurrence(T,sp,kr .- 1,tocanonical(sp,x))
332+
else
333+
z = Zeros{T}(length(range(minimum(kr), order, step=step(kr))))
334+
kr_red = kr .- (order + 1)
335+
polydegrees = reverse(range(maximum(kr_red), max(0, minimum(kr_red)), step=-step(kr)))
336+
D = Derivative(sp, order)
337+
if !isempty(polydegrees)
338+
P = forwardrecurrence(T, rangespace(D), polydegrees, tocanonical(sp,x))
339+
# in case the derivative has only one non-zero band, the matrix-vector
340+
# multiplication simplifies considerably.
341+
# This branch is particularly useful for ConcreteDerivatives where
342+
# indexing is fast, as the non-zero band may be computed without
343+
# evaluating the matrix
344+
if bandwidth(D, 1) == -bandwidth(D, 2)
345+
d = nonzeroband(T, D, polydegrees, order)
346+
Pd = P .* d
347+
if !isempty(z)
348+
Pd = T[z; Pd]
349+
end
350+
else
351+
# in general, this is a banded matrix-vector product
352+
Dtv = view(transpose(D), kr, 1:length(P))
353+
Pd = strictconvert(Vector{T}, mul_coefficients(Dtv, P))
354+
end
355+
else
356+
Pd = T[z;]
357+
end
358+
Pd
359+
end
360+
end
361+
362+
function nonzeroband(::Type{T}, D::ConcreteDerivative, polydegrees, order) where {T}
363+
bw = bandwidth(D, 2)
364+
T[D[k+1, k+1+bw] for k in polydegrees]
365+
end
366+
function nonzeroband(::Type{T}, D, polydegrees, order) where {T}
367+
bw = bandwidth(D, 2)
368+
rows = 1:(maximum(polydegrees) + order + 1)
369+
B = D[rows, rows .+ bw]
370+
Bv = @view B[diagind(B)]
371+
Bv[polydegrees .+ 1]
337372
end
338373

339374

0 commit comments

Comments
 (0)