Skip to content

Commit 943b781

Browse files
dlfivefiftydlfivefifty
committed
Revert "try using chebyshevtransform instead of applyTN"
This reverts commit 02cba99.
1 parent 02cba99 commit 943b781

11 files changed

+23
-66
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₂ = plan_ichebyshevtransform!(c₁; kind=2), plan_DSTI!(c₂)
130+
p₁,p₂ = applyTN_plan(c₁),applyUN_plan(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₂ = plan_ichebyshevtransform!(c₁; kind=2),plan_DSTI!(c₂)
171+
p₁,p₂ = applyTN_plan(c₁),applyUN_plan(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₂ = plan_ichebyshevtransform!(c₁; kind=2),plan_DSTI!(c₂)
128+
p₁,p₂ = applyTN_plan(c₁),applyUN_plan(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₂ = plan_ichebyshevtransform!(c₁; kind=2),plan_DSTI!(c₂)
168+
p₁,p₂ = applyTN_plan(c₁),applyUN_plan(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: 3 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,7 @@ 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, plan_r2r!,
11-
REDFT00, REDFT01, REDFT10, REDFT11,
12-
RODFT00, RODFT01, RODFT10, RODFT11
10+
import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan
1311
const LAmul! = Base.A_mul_B!
1412
import Base: Factorization
1513
rmul!(A::AbstractArray, c::Number) = scale!(A,c)
@@ -19,9 +17,7 @@ if VERSION < v"0.7-"
1917
else
2018
using FFTW, LinearAlgebra, DSP
2119
import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
22-
import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan, plan_r2r!,
23-
REDFT00, REDFT01, REDFT10, REDFT11,
24-
RODFT00, RODFT01, RODFT10, RODFT11
20+
import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan
2521
const LAmul! = LinearAlgebra.mul!
2622
const libfftw = libfftw3
2723
const libfftwf = libfftw3f
@@ -30,7 +26,7 @@ else
3026
end
3127

3228

33-
import Base: *, \, inv, size, view
29+
import Base: *, \, size, view
3430
import Base: getindex, setindex!, length
3531
import Compat.LinearAlgebra: BlasFloat, BlasInt
3632
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-
p₁*c_cheb2
14+
applyTN!(c_cheb2,p₁)
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-
p₁*c₁; p₂*c₂
46+
applyTN!(c₁,p₁);applyUN!(c₂,p₂)
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-
p₁*c_cheb2
14+
applyTN!(c_cheb2,p₁)
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-
p₁*c₁; p₂*c₂
47+
applyTN!(c₁,p₁);applyUN!(c₂,p₂)
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: 5 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,5 @@
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
1+
export plan_chebyshevtransform, plan_ichebyshevtransform, chebyshevtransform, ichebyshevtransform, chebyshevpoints,
2+
plan_chebyshevutransform, plan_ichebyshevutransform, chebyshevutransform, ichebyshevutransform
53

64
## Transforms take values at Chebyshev points of the first and second kinds and produce Chebyshev coefficients
75

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

6664
chebyshevtransform(x;kind::Integer=1) = chebyshevtransform!(copy(x);kind=kind)
6765

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

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

@@ -77,15 +72,6 @@ struct IChebyshevTransformPlan{T,kind,inplace,P}
7772
plan::P
7873
end
7974

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-
8975
function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
9076
if kind == 1
9177
if length(x) == 0
@@ -97,7 +83,8 @@ function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T
9783
if length(x) ≤ 1
9884
error("Cannot create a length $(length(x)) inverse chebyshev transform")
9985
end
100-
inv(plan_chebyshevtransform!(x;kind=2))
86+
plan = plan_chebyshevtransform!(x;kind=2)
87+
IChebyshevTransformPlan{T,2,true,typeof(plan)}(plan)
10188
end
10289
end
10390

@@ -320,22 +307,3 @@ function chebyshevpoints(::Type{T}, n::Integer; kind::Int=1) where T<:Number
320307
end
321308
end
322309
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-
p₁*c₁; p₂*c₂
23+
applyTN!(c₁,p₁);applyUN!(c₂,p₂)
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-
p₁ \ v_cheb
40+
applyTNinv!(v_cheb,p₁)
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-
p₁*c₁; p₂*c₂
24+
applyTN!(c₁,p₁);applyUN!(c₂,p₂)
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-
p₁ \ v_cheb
45+
applyTNinv!(v_cheb,p₁)
4646
end

‎test/chebyshevtests.jl

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -23,11 +23,6 @@ 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))
3126
end
3227
end
3328

‎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)
96
end
107

118
@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)