Inplace jacobi/Legendre transforms (#81)

* inplace jacobi/Legendre transforms

* tests for non-zero orders

* tests for fractional orders

* move inplace flag to type-domain

Author: jishnub

src/Spaces/Jacobi/jacobitransform.jl

@@ -11,22 +11,31 @@ JacobiTransformPlan(chebplan::CPLAN, cjtplan::CJT) where {CPLAN,CJT} =
 
 plan_transform(S::Jacobi, v::AbstractVector) =
     JacobiTransformPlan(plan_transform(Chebyshev(), v), plan_cheb2jac(v, S.a, S.b))
-*(P::JacobiTransformPlan, vals::AbstractVector) = P.cjtplan*(P.chebplan*vals)
+plan_transform!(S::Jacobi, v::AbstractVector) =
+    JacobiTransformPlan(plan_transform!(Chebyshev(), v), plan_cheb2jac(v, S.a, S.b))
+*(P::JacobiTransformPlan, vals::AbstractVector) = lmul!(P.cjtplan, P.chebplan * vals)
 
 
-struct JacobiITransformPlan{T,CPLAN,CJT} <: AbstractTransformPlan{T}
+struct JacobiITransformPlan{T,CPLAN,CJT,inplace} <: AbstractTransformPlan{T}
     ichebplan::CPLAN
     icjtplan::CJT
 end
 
-JacobiITransformPlan(chebplan::CPLAN, cjtplan::CJT) where {CPLAN,CJT} =
-    JacobiITransformPlan{eltype(chebplan),CPLAN,CJT}(chebplan, cjtplan)
-
+JacobiITransformPlan(chebplan, cjtplan) = JacobiITransformPlan{false}(chebplan, cjtplan)
+JacobiITransformPlan{inplace}(chebplan::CPLAN, cjtplan::CJT) where {CPLAN,CJT,inplace} =
+    JacobiITransformPlan{eltype(chebplan),CPLAN,CJT,inplace}(chebplan, cjtplan)
 
+inplace(J::JacobiITransformPlan{<:Any,<:Any,<:Any,IP}) where {IP} = IP
 
 plan_itransform(S::Jacobi, v::AbstractVector) =
     JacobiITransformPlan(plan_itransform(Chebyshev(), v), plan_jac2cheb(v, S.a, S.b))
-*(P::JacobiITransformPlan, cfs::AbstractVector) = P.ichebplan*(P.icjtplan*cfs)
+plan_itransform!(S::Jacobi, v::AbstractVector) =
+    JacobiITransformPlan{true}(plan_itransform!(Chebyshev(), v), plan_jac2cheb(v, S.a, S.b))
+icjt(P, cfs, ::Val{true}) = lmul!(P, cfs)
+icjt(P, cfs, ::Val{false}) = P * cfs
+function *(P::JacobiITransformPlan, cfs::AbstractVector)
+    P.ichebplan * icjt(P.icjtplan, cfs, Val(inplace(P)))
+end
 
 
 function coefficients(f::AbstractVector{T}, a::Jacobi, b::Chebyshev) where T

test/JacobiTest.jl

@@ -1,7 +1,7 @@
 using ApproxFunOrthogonalPolynomials, ApproxFunBase, Test, SpecialFunctions, LinearAlgebra
 import ApproxFunBase: testbandedbelowoperator, testbandedoperator, testspace, testtransforms, Vec,
                     maxspace, NoSpace, hasconversion, testfunctional,
-                    reverseorientation, ReverseOrientation
+                    reverseorientation, ReverseOrientation, transform!, itransform!
 import ApproxFunOrthogonalPolynomials: jacobip
 
 @testset "Jacobi" begin
@@ -33,6 +33,21 @@ import ApproxFunOrthogonalPolynomials: jacobip
     @testset "Conversion" begin
         testtransforms(Jacobi(-0.5,-0.5))
         @test norm(Fun(Fun(exp),Jacobi(-.5,-.5))-Fun(exp,Jacobi(-.5,-.5))) < 300eps()
+
+        @testset "inplace transform" begin
+            for T in [Float32, Float64, BigFloat]
+                for v in Any[rand(T, 10), rand(complex(T), 10)]
+                    v2 = copy(v)
+                    for a in 0:0.5:3, b in 0:0.5:3
+                        J = Jacobi(a, b)
+                        transform!(J, v)
+                        @test transform(J, v2) ≈ v
+                        itransform!(J, v)
+                        @test v2 ≈ v
+                    end
+                end
+            end
+        end
     end
 
     @testset "Derivative" begin