Improvements in Legendre Transforms (#124)

* Work on matrix Legendre transforms

* LegendreTransformPlan

* tests pass

* v0.7.3

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.7.2"
+version = "0.7.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -33,7 +33,7 @@ ContinuumArrays = "0.12.4"
 DomainSets = "0.5.6, 0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3, 0.5"
-FastTransforms = "0.14.9"
+FastTransforms = "0.15"
 FillArrays = "0.13"
 HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.11"

src/classical/legendre.jl

@@ -79,11 +79,28 @@ end
 
 
 ldiv(P::Legendre{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)
-function ldiv(::Legendre{V}, f::AbstractQuasiVector) where V
+function transform_ldiv(::Legendre{V}, f::Union{AbstractQuasiVector,AbstractQuasiMatrix}) where V
     T = ChebyshevT{V}()
-    [cheb2leg(paddeddata(T \ f)); zeros(V,∞)]
+    dat = T \ f
+    pad(cheb2leg(paddeddata(dat)), axes(dat)...)
 end
 
+struct LegendreTransformPlan{T, CHEB2LEG, DCT} <: Plan{T}
+    cheb2leg::CHEB2LEG
+    chebtransform::DCT
+end
+
+LegendreTransformPlan(c2l, ct) = LegendreTransformPlan{promote_type(eltype(c2l),eltype(ct)),typeof(c2l),typeof(ct)}(c2l, ct)
+
+*(P::LegendreTransformPlan, x::AbstractArray) = P.cheb2leg*(P.chebtransform*x)
+
+function plan_grid_transform(P::Legendre{T}, szs::NTuple{N,Int}, dims=1:N) where {T,N}
+    arr = Array{T}(undef, szs...)
+    x = grid(P, size(arr,1))
+    x, LegendreTransformPlan(FastTransforms.plan_th_cheb2leg!(arr, dims), plan_chebyshevtransform(arr, dims))
+end
+
+
 """
     legendre_massmatrix
 

test/test_chebyshev.jl

@@ -121,6 +121,18 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
             @test U[0.1,:] ≈ cos.(0.1 * (0:1000))
             @test U[[0.1,0.2],:] ≈ [cos.(0.1 * (0:1000)'); cos.(0.2 * (0:1000)')]
 
+            # plan_transform for rect
+
+            X = randn(10, 11)
+            F1 = plan_transform(T, X, 1)
+            p1 = plan_transform(T, X[:,1])
+            @test F1*X ≈ hcat([p1 * X[:,j] for j = 1:size(X,2)]...)
+            F2 = plan_transform(T, X, 2)
+            p2 = plan_transform(T, X[1,:])
+            @test F2*X ≈ vcat([(p2 * X[k,:])' for k = 1:size(X,1)]...)
+            F = plan_transform(T, X)
+            @test F*X ≈ F2*(F1*X)
+
             # support tensors but for grids
             X = randn(150, 2, 2)
             F = plan_transform(T, X, 1)

test/test_legendre.jl

@@ -40,6 +40,24 @@ import QuasiArrays: MulQuasiArray
         end
     end
 
+    @testset "expansion" begin
+        P = Legendre()
+        x = axes(P,1)
+        @test (P \ [exp.(x) cos.(x)])[1:10,1:2] ≈ [P\exp.(x) P\cos.(x)][1:10,:]
+
+        # plan_transform for rect
+
+        X = randn(10, 11)
+        F1 = plan_transform(P, X, 1)
+        p1 = plan_transform(P, X[:,1])
+        @test F1*X ≈ hcat([p1 * X[:,j] for j = 1:size(X,2)]...)
+        F2 = plan_transform(P, X, 2)
+        p2 = plan_transform(P, X[1,:])
+        @test F2*X ≈ vcat([(p2 * X[k,:])' for k = 1:size(X,1)]...)
+        F = plan_transform(P, X)
+        @test F*X ≈ F2*(F1*X)
+    end
+
     @testset "operators" begin
         P = Legendre()
         P̃ = Jacobi(0.0,0.0)