Skip to content

Commit 02cba99

Browse files
dlfivefiftydlfivefifty
committed
try using chebyshevtransform instead of applyTN
1 parent 08a429b commit 02cba99

11 files changed

+66
-23
lines changed

‎src/ChebyshevJacobiPlan.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -127,7 +127,7 @@ function ForwardChebyshevJacobiPlan(c_jac::AbstractVector{T},α::T,β::T,M::Int)
127127
c₁,c₂,um,vm = zero(c_jac),view(zero(c_jac),2:N),zero(c_jac),zero(c_jac)
128128

129129
# Initialize DCT-I and DST-I plans
130-
p₁,p₂ = applyTN_plan(c₁),applyUN_plan(c₂)
130+
p₁,p₂ = plan_ichebyshevtransform!(c₁; kind=2), plan_DSTI!(c₂)
131131

132132
# Initialize coefficients of the asymptotic formula
133133
cfs = init_cfs(α,β,M)
@@ -168,7 +168,7 @@ function BackwardChebyshevJacobiPlan(c_cheb::AbstractVector{T},α::T,β::T,M::In
168168
c₁,c₂,um,vm = zero(c_cheb2),view(zero(c_cheb2),2:2N),zero(c_cheb2),zero(c_cheb2)
169169

170170
# Initialize DCT-I and DST-I plans
171-
p₁,p₂ = applyTN_plan(c₁),applyUN_plan(c₂)
171+
p₁,p₂ = plan_ichebyshevtransform!(c₁; kind=2),plan_DSTI!(c₂)
172172

173173
# Initialize coefficients of the asymptotic formula
174174
cfs = init_cfs(α,β,M)

‎src/ChebyshevUltrasphericalPlan.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -125,7 +125,7 @@ function ForwardChebyshevUltrasphericalPlan(c_ultra::AbstractVector{T},λ::T,M::
125125
c₁,c₂,um,vm = zero(c_ultra),view(zero(c_ultra),2:N),zero(c_ultra),zero(c_ultra)
126126

127127
# Initialize DCT-I and DST-I plans
128-
p₁,p₂ = applyTN_plan(c₁),applyUN_plan(c₂)
128+
p₁,p₂ = plan_ichebyshevtransform!(c₁; kind=2),plan_DSTI!(c₂)
129129

130130
# Clenshaw-Curtis points
131131
θ = N > 0 ? T[k/N for k=zero(T):N] : T[0]
@@ -165,7 +165,7 @@ function BackwardChebyshevUltrasphericalPlan(c_ultra::AbstractVector{T},λ::T,M:
165165
c₁,c₂,um,vm = zero(c_cheb2),view(zero(c_cheb2),2:2N),zero(c_cheb2),zero(c_cheb2)
166166

167167
# Initialize DCT-I and DST-I plans
168-
p₁,p₂ = applyTN_plan(c₁),applyUN_plan(c₂)
168+
p₁,p₂ = plan_ichebyshevtransform!(c₁; kind=2),plan_DSTI!(c₂)
169169

170170
# Clenshaw-Curtis nodes and weights
171171
θ = N > 0 ? T[k/2N for k=zero(T):2N] : T[0]

‎src/FastTransforms.jl

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,9 @@ using ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter, Comp
77
if VERSION < v"0.7-"
88
using Base.FFTW
99
import Base.FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
10-
import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan
10+
import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan, plan_r2r!,
11+
REDFT00, REDFT01, REDFT10, REDFT11,
12+
RODFT00, RODFT01, RODFT10, RODFT11
1113
const LAmul! = Base.A_mul_B!
1214
import Base: Factorization
1315
rmul!(A::AbstractArray, c::Number) = scale!(A,c)
@@ -17,7 +19,9 @@ if VERSION < v"0.7-"
1719
else
1820
using FFTW, LinearAlgebra, DSP
1921
import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
20-
import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan
22+
import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan, plan_r2r!,
23+
REDFT00, REDFT01, REDFT10, REDFT11,
24+
RODFT00, RODFT01, RODFT10, RODFT11
2125
const LAmul! = LinearAlgebra.mul!
2226
const libfftw = libfftw3
2327
const libfftwf = libfftw3f
@@ -26,7 +30,7 @@ else
2630
end
2731

2832

29-
import Base: *, \, size, view
33+
import Base: *, \, inv, size, view
3034
import Base: getindex, setindex!, length
3135
import Compat.LinearAlgebra: BlasFloat, BlasInt
3236
import HierarchicalMatrices: HierarchicalMatrix, unsafe_broadcasttimes!

‎src/cheb2jac.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ function cheb2jac(c_cheb::AbstractVector{T},α::T,β::T,plan::ChebyshevJacobiPla
1111
v_jac = zero(c_cheb2)
1212

1313
# Perform transposed DCT-I
14-
applyTN!(c_cheb2,p₁)
14+
p₁*c_cheb2
1515

1616
# Scale values by Clenshaw-Curtis weights
1717
@inbounds for i=1:2N+1 c_cheb2[i] *= w[i] end
@@ -43,7 +43,7 @@ function cheb2jac(c_cheb::AbstractVector{T},α::T,β::T,plan::ChebyshevJacobiPla
4343
init_c₁c₂!(c₁,c₂.parent,um,vm,c_cheb2,i₁[k+1],i₂[k+1])
4444

4545
# Apply planned DCT-I and DST-I in-place
46-
applyTN!(c₁,p₁);applyUN!(c₂,p₂)
46+
p₁*c₁; p₂*c₂
4747

4848
# Compute diagonal 2N-scaling multiplied by local coefficients and zero out excess
4949
@inbounds for j=j₁[k]:j₂[k] v_jac[j] += cnmαβ[j]*(c₁[j]+c₂.parent[j]) end

‎src/cheb2ultra.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ function cheb2ultra(c_cheb::AbstractVector{T},λ::T,plan::ChebyshevUltraspherica
1111
v_ultra = zero(c_cheb2)
1212

1313
# Perform transposed DCT-I
14-
applyTN!(c_cheb2,p₁)
14+
p₁*c_cheb2
1515

1616
# Scale values by Clenshaw-Curtis weights
1717
@inbounds for i=1:2N+1 c_cheb2[i] *= w[i] end
@@ -44,7 +44,7 @@ function cheb2ultra(c_cheb::AbstractVector{T},λ::T,plan::ChebyshevUltraspherica
4444
init_c₁c₂!(c₁,c₂.parent,um,vm,c_cheb2,i₁[k+1],i₂[k+1])
4545

4646
# Apply planned DCT-I and DST-I in-place
47-
applyTN!(c₁,p₁);applyUN!(c₂,p₂)
47+
p₁*c₁; p₂*c₂
4848

4949
# Compute diagonal 2N-scaling multiplied by local coefficients and zero out excess
5050
@inbounds for j=j₁[k]:j₂[k] v_ultra[j] += cnmλ[j]*(c₁[j]+c₂.parent[j]) end

‎src/chebyshevtransform.jl

Lines changed: 37 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,7 @@
1-
export plan_chebyshevtransform, plan_ichebyshevtransform, chebyshevtransform, ichebyshevtransform, chebyshevpoints,
2-
plan_chebyshevutransform, plan_ichebyshevutransform, chebyshevutransform, ichebyshevutransform
1+
export plan_chebyshevtransform, plan_ichebyshevtransform, plan_chebyshevtransform!, plan_ichebyshevtransform!,
2+
chebyshevtransform, ichebyshevtransform, chebyshevpoints,
3+
plan_chebyshevutransform, plan_ichebyshevutransform, plan_chebyshevutransform!, plan_ichebyshevutransform!,
4+
chebyshevutransform, ichebyshevutransform
35

46
## Transforms take values at Chebyshev points of the first and second kinds and produce Chebyshev coefficients
57

@@ -63,7 +65,10 @@ chebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber}
6365

6466
chebyshevtransform(x;kind::Integer=1) = chebyshevtransform!(copy(x);kind=kind)
6567

66-
*(P::ChebyshevTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
68+
*(P::ChebyshevTransformPlan{T,k,false}, x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
69+
70+
71+
6772

6873
## Inverse transforms take Chebyshev coefficients and produce values at Chebyshev points of the first and second kinds
6974

@@ -72,6 +77,15 @@ struct IChebyshevTransformPlan{T,kind,inplace,P}
7277
plan::P
7378
end
7479

80+
# second kind Chebyshev transforms share a plan with their inverse
81+
# so we support this via inv
82+
inv(P::ChebyshevTransformPlan{T,2,true}) where T = IChebyshevTransformPlan{T,2,true,typeof(P)}(P)
83+
inv(P::IChebyshevTransformPlan{T,2,true}) where T = P.plan
84+
85+
\(P::ChebyshevTransformPlan, x::AbstractArray) = inv(P) * x
86+
\(P::IChebyshevTransformPlan, x::AbstractArray) = inv(P) * x
87+
88+
7589
function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
7690
if kind == 1
7791
if length(x) == 0
@@ -83,8 +97,7 @@ function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T
8397
if length(x) ≤ 1
8498
error("Cannot create a length $(length(x)) inverse chebyshev transform")
8599
end
86-
plan = plan_chebyshevtransform!(x;kind=2)
87-
IChebyshevTransformPlan{T,2,true,typeof(plan)}(plan)
100+
inv(plan_chebyshevtransform!(x;kind=2))
88101
end
89102
end
90103

@@ -307,3 +320,22 @@ function chebyshevpoints(::Type{T}, n::Integer; kind::Int=1) where T<:Number
307320
end
308321
end
309322
chebyshevpoints(n::Integer; kind::Int=1) = chebyshevpoints(Float64, n; kind=kind)
323+
324+
325+
# sin(nθ) coefficients to values at Clenshaw-Curtis nodes except ±1
326+
327+
struct DSTPlan{T,kind,inplace,P} <: Plan{T}
328+
plan::P
329+
end
330+
331+
DSTPlan{k,inp}(plan) where {k,inp} =
332+
DSTPlan{eltype(plan),k,inp,typeof(plan)}(plan)
333+
334+
335+
plan_DSTI!(x) = length(x) > 0 ? DSTPlan{1,true}(FFTW.plan_r2r!(x, FFTW.RODFT00)) :
336+
ones(x)'
337+
338+
function *(P::DSTPlan{T,1}, x::AbstractArray) where {T}
339+
x = P.plan*x
340+
rmul!(x,half(T))
341+
end

‎src/jac2cheb.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ function jac2cheb(c_jac::AbstractVector{T},α::T,β::T,plan::ChebyshevJacobiPlan
2020
init_c₁c₂!(c₁,c₂.parent,cnmαβ,c_jac,j₁[k],j₂[k])
2121

2222
# Apply planned DCT-I and DST-I in-place
23-
applyTN!(c₁,p₁);applyUN!(c₂,p₂)
23+
p₁*c₁; p₂*c₂
2424

2525
# Compute u_m(θ) and v_m(θ)
2626
compute_umvm!(um,vm,cfs,α,β,tempcos,tempsin,tempcosβsinα,m,θ,i₁[k+1]:i₂[k+1])
@@ -37,5 +37,5 @@ function jac2cheb(c_jac::AbstractVector{T},α::T,β::T,plan::ChebyshevJacobiPlan
3737
@inbounds for i=i₁[k+1]:i₂[k+1] v_cheb[i] += backward_recurrence(c_jac,j₂[k+1],θ[i],tempcos[i],tempsin[i],rp) end
3838

3939
# perform IDCT-I
40-
applyTNinv!(v_cheb,p₁)
40+
p₁ \ v_cheb
4141
end

‎src/ultra2cheb.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ function ultra2cheb(c_ultra::AbstractVector{T},λ::T,plan::ChebyshevUltraspheric
2121
init_c₁c₂!(c₁,c₂.parent,cnmλ,c_ultra,j₁[k],j₂[k])
2222

2323
# Apply planned DCT-I and DST-I in-place
24-
applyTN!(c₁,p₁);applyUN!(c₂,p₂)
24+
p₁*c₁; p₂*c₂
2525

2626
# Compute u_m(θ) and v_m(θ)
2727
compute_umvm!(um,vm,λ,tempsinλm,m,θ,i₁[k+1]:i₂[k+1])
@@ -42,5 +42,5 @@ function ultra2cheb(c_ultra::AbstractVector{T},λ::T,plan::ChebyshevUltraspheric
4242
@inbounds for i=i₁[k+1]:i₂[k+1] v_cheb[i] += backward_recurrence(c_ultra,j₂[k+1],θ[i],tempsin2[N+2-i],tempsin2[i],rp) end
4343

4444
# perform IDCT-I
45-
applyTNinv!(v_cheb,p₁)
45+
p₁ \ v_cheb
4646
end

‎test/chebyshevtests.jl

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,11 @@ using FastTransforms, Compat, Compat.Test
2323
f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
2424
@test f̃(0.1) ≈ exp(T(0.1))
2525
@test ichebyshevtransform(f̌; kind=2) ≈ exp.(p_2)
26+
27+
P = plan_chebyshevtransform!(f; kind=2)
28+
Pi = plan_ichebyshevtransform!(f; kind=2)
29+
@test all((P \ copy(f)) .=== Pi * copy(f))
30+
@test f ≈ P \ (P*copy(f)) ≈ P * (P\copy(f))
2631
end
2732
end
2833

‎test/nuffttests.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3,12 +3,12 @@ using Compat.Test
33

44
if VERSION ≥ v"0.7-"
55
using FFTW
6+
FFTW.set_num_threads(Base.Sys.CPU_THREADS)
7+
else
8+
FFTW.set_num_threads(Base.Sys.CPU_CORES)
69
end
710

811
@testset "Nonuniform fast Fourier transforms" begin
9-
10-
FFTW.set_num_threads(Base.Sys.CPU_CORES)
11-
1212
function nudft1(c::AbstractVector, ω::AbstractVector{T}) where {T<:AbstractFloat}
1313
# Nonuniform discrete Fourier transform of type I
1414

0 commit comments

Comments
 (0)