Just use applyTN for now

Author: dlfivefifty

src/FastTransforms.jl

@@ -7,7 +7,9 @@ using ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter, Comp
 if VERSION < v"0.7-"
     using Base.FFTW
     import Base.FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
-    import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan
+    import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan, plan_r2r!,
+                        REDFT00, REDFT01, REDFT10, REDFT11,
+                        RODFT00, RODFT01, RODFT10, RODFT11
     const LAmul! = Base.A_mul_B!
     import Base: Factorization
     rmul!(A::AbstractArray, c::Number) = scale!(A,c)
@@ -17,7 +19,9 @@ if VERSION < v"0.7-"
 else
     using FFTW, LinearAlgebra, DSP
     import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
-    import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan
+    import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan, plan_r2r!,
+                    REDFT00, REDFT01, REDFT10, REDFT11,
+                    RODFT00, RODFT01, RODFT10, RODFT11
     const LAmul! = LinearAlgebra.mul!
     const libfftw = libfftw3
     const libfftwf = libfftw3f
@@ -26,7 +30,7 @@ else
 end
 
 
-import Base: *, \, size, view
+import Base: *, \, inv, size, view
 import Base: getindex, setindex!, length
 import Compat.LinearAlgebra: BlasFloat, BlasInt
 import HierarchicalMatrices: HierarchicalMatrix, unsafe_broadcasttimes!

src/chebyshevtransform.jl

@@ -1,5 +1,7 @@
-export plan_chebyshevtransform, plan_ichebyshevtransform, chebyshevtransform, ichebyshevtransform, chebyshevpoints,
-            plan_chebyshevutransform, plan_ichebyshevutransform, chebyshevutransform, ichebyshevutransform
+export plan_chebyshevtransform, plan_ichebyshevtransform, plan_chebyshevtransform!, plan_ichebyshevtransform!,
+            chebyshevtransform, ichebyshevtransform, chebyshevpoints,
+            plan_chebyshevutransform, plan_ichebyshevutransform, plan_chebyshevutransform!, plan_ichebyshevutransform!,
+            chebyshevutransform, ichebyshevutransform
 
 ## Transforms take values at Chebyshev points of the first and second kinds and produce Chebyshev coefficients
 
@@ -63,7 +65,10 @@ chebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber}
 
 chebyshevtransform(x;kind::Integer=1) = chebyshevtransform!(copy(x);kind=kind)
 
-*(P::ChebyshevTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
+*(P::ChebyshevTransformPlan{T,k,false}, x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
+
+
+
 
 ## Inverse transforms take Chebyshev coefficients and produce values at Chebyshev points of the first and second kinds
 
@@ -72,6 +77,15 @@ struct IChebyshevTransformPlan{T,kind,inplace,P}
     plan::P
 end
 
+# second kind Chebyshev transforms share a plan with their inverse
+# so we support this via inv
+inv(P::ChebyshevTransformPlan{T,2,true}) where T = IChebyshevTransformPlan{T,2,true,typeof(P)}(P)
+inv(P::IChebyshevTransformPlan{T,2,true}) where T = P.plan
+
+\(P::ChebyshevTransformPlan, x::AbstractArray) = inv(P) * x
+\(P::IChebyshevTransformPlan, x::AbstractArray) = inv(P) * x
+
+
 function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
     if kind == 1
         if length(x) == 0
@@ -83,8 +97,7 @@ function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T
         if length(x) ≤ 1
             error("Cannot create a length $(length(x)) inverse chebyshev transform")
         end
-        plan = plan_chebyshevtransform!(x;kind=2)
-        IChebyshevTransformPlan{T,2,true,typeof(plan)}(plan)
+        inv(plan_chebyshevtransform!(x;kind=2))
     end
 end
 
@@ -307,3 +320,22 @@ function chebyshevpoints(::Type{T}, n::Integer; kind::Int=1) where T<:Number
     end
 end
 chebyshevpoints(n::Integer; kind::Int=1) = chebyshevpoints(Float64, n; kind=kind)
+
+
+# sin(nθ) coefficients to values at Clenshaw-Curtis nodes except ±1
+
+struct DSTPlan{T,kind,inplace,P} <: Plan{T}
+    plan::P
+end
+
+DSTPlan{k,inp}(plan) where {k,inp} =
+    DSTPlan{eltype(plan),k,inp,typeof(plan)}(plan)
+
+
+plan_DSTI!(x) = length(x) > 0 ? DSTPlan{1,true}(FFTW.plan_r2r!(x, FFTW.RODFT00)) :
+                                ones(x)'
+
+function *(P::DSTPlan{T,1}, x::AbstractArray) where {T}
+    x = P.plan*x
+    rmul!(x,half(T))
+end

src/clenshawcurtis.jl

@@ -17,3 +17,46 @@ function clenshawcurtisweights!(μ::Vector{T}, plan) where T
     μ[1] *= half(T); μ[N] *= half(T)
     return μ
 end
+
+
+
+# Chebyshev-T coefficients to values at Clenshaw-Curtis nodes
+
+applyTN_plan(x) = length(x) > 1 ? FFTW.plan_r2r!(x, FFTW.REDFT00) : fill!(similar(x),1)'
+
+applyTN!(x::Vector{T}) where {T<:AbstractFloat} = applyTN!(x,applyTN_plan(x))
+
+function applyTN!(x::Vector{T},plan) where T<:AbstractFloat
+    x[1] *= 2; x[end] *=2
+    plan*x
+    rmul!(x,half(T))
+end
+applyTN(x::Vector{T},plan) where {T<:AbstractFloat} = applyTN!(copy(x),plan)
+applyTN(x::Vector{T}) where {T<:AbstractFloat} = applyTN!(copy(x))
+
+# Values at Clenshaw-Curtis nodes to Chebyshev-T coefficients
+
+applyTNinv_plan(x) = length(x) > 1 ? FFTW.plan_r2r!(x, FFTW.REDFT00) : ones(x)'
+
+applyTNinv!(x::Vector{T}) where {T<:AbstractFloat} = applyTNinv!(x,applyTNinv_plan(x))
+
+function applyTNinv!(x::Vector{T},plan) where T<:AbstractFloat
+    plan*x
+    x[1] /= 2;x[end] /= 2
+    rmul!(x,inv(length(x)-one(T)))
+end
+applyTNinv(x::Vector{T},plan) where {T<:AbstractFloat} = applyTNinv!(copy(x),plan)
+applyTNinv(x::Vector{T}) where {T<:AbstractFloat} = applyTNinv!(copy(x))
+
+# sin(nθ) coefficients to values at Clenshaw-Curtis nodes except ±1
+
+applyUN_plan(x) = length(x) > 0 ? FFTW.plan_r2r!(x, FFTW.RODFT00) : fill!(similar(x),1)'
+
+applyUN!(x::AbstractVector{T}) where {T<:AbstractFloat} = applyUN!(x,applyUN_plan(x))
+
+function applyUN!(x::AbstractVector{T},plan) where T<:AbstractFloat
+    plan*x
+    rmul!(x,half(T))
+end
+applyUN(x::AbstractVector{T},plan) where {T<:AbstractFloat} = applyUN!(copy(x),plan)
+applyUN(x::AbstractVector{T}) where {T<:AbstractFloat} = applyUN!(copy(x))

test/chebyshevtests.jl

@@ -23,6 +23,11 @@ using FastTransforms, Compat, Compat.Test
             f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
             @test f̃(0.1) ≈ exp(T(0.1))
             @test ichebyshevtransform(f̌; kind=2) ≈ exp.(p_2)
+
+            P = plan_chebyshevtransform!(f; kind=2)
+            Pi = plan_ichebyshevtransform!(f; kind=2)
+            @test all((P \ copy(f)) .=== Pi * copy(f))
+            @test f ≈ P \ (P*copy(f)) ≈ P * (P\copy(f))
         end
     end
 

test/nuffttests.jl

@@ -3,12 +3,12 @@ using Compat.Test
 
 if VERSION ≥ v"0.7-"
     using FFTW
+    FFTW.set_num_threads(Base.Sys.CPU_THREADS)
+else
+    FFTW.set_num_threads(Base.Sys.CPU_CORES)
 end
 
 @testset "Nonuniform fast Fourier transforms" begin
-
-    FFTW.set_num_threads(Base.Sys.CPU_CORES)
-
     function nudft1(c::AbstractVector, ω::AbstractVector{T}) where {T<:AbstractFloat}
         # Nonuniform discrete Fourier transform of type I
 

test/runtests.jl

@@ -5,6 +5,8 @@ srand(0)
 
 include("basictests.jl")
 
+include("chebyshevtests.jl")
+
 include("nuffttests.jl")
 
 include("chebyshevjacobitests.jl")

test/toeplitztests.jl

@@ -1,6 +1,6 @@
 using FastTransforms, ToeplitzMatrices
 
-@testset "BigFloat TOeplitz" begin
+@testset "BigFloat Toeplitz" begin
     T = Toeplitz(BigFloat[1,2,3,4,5], BigFloat[1,6,7,8,0])
     @test T*ones(BigFloat,5) ≈ [22,24,19,16,15]