Support 2D Clenshaw on rectangles

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.7.2"
+version = "0.7.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -19,6 +19,8 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+RecurrenceRelationshipArrays = "b889d2dc-af3c-4820-88a8-238fa91d3518"
+RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -31,7 +33,7 @@ BlockBandedMatrices = "0.13"
 ClassicalOrthogonalPolynomials = "0.14"
 ContinuumArrays = "0.18"
 DomainSets = "0.6, 0.7"
-FastTransforms = "0.16"
+FastTransforms = "0.17"
 FillArrays = "1.0"
 HarmonicOrthogonalPolynomials = "0.6"
 InfiniteArrays = "0.15"

src/rect.jl

@@ -148,7 +148,7 @@ function Base.unsafe_getindex(f::Mul{MultivariateOPLayout{2},<:DiagTravLayout{<:
     P,c = f.A, f.B
     A,B = P.args
     x,y = 𝐱
-    clenshaw(clenshaw(paddeddata(c.array), recurrencecoefficients(A)..., x), recurrencecoefficients(B)..., y)
+    clenshaw(vec(clenshaw(paddeddata(c.array), recurrencecoefficients(A)..., x; dims=1)), recurrencecoefficients(B)..., y)
 end
 
 Base.@propagate_inbounds function getindex(f::Mul{MultivariateOPLayout{2},<:DiagTravLayout{<:PaddedLayout}}, x::SVector, j...)

test/test_rect.jl

@@ -1,5 +1,5 @@
 using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, Test
-using ClassicalOrthogonalPolynomials: expand
+using ClassicalOrthogonalPolynomials: expand, coefficients, recurrencecoefficients
 using MultivariateOrthogonalPolynomials: weaklaplacian
 using ContinuumArrays: plotgridvalues
 
@@ -154,4 +154,33 @@ using ContinuumArrays: plotgridvalues
         Dₓ = KronTrav(D²,M)
         @test Dₓ[Block.(1:1),Block.(1:1)] == Dₓ[Block(1,1)]
     end
+
+    @testset "variable coefficients" begin
+        T,U = ChebyshevT(), ChebyshevU()
+        P = RectPolynomial(T, U)
+        𝐱 = axes(P,1)
+        x,y = first.(𝐱), last.(𝐱)
+        X = P\(x .* P)
+        Y = P\(y .* P)
+
+        @test X isa KronTrav
+        @test Y isa KronTrav
+
+        a =  (x,y) -> I + x + 2y + 3x^2 +4x*y + 5y^2
+        𝐚 = expand(P,splat(a))
+        A = a(X,Y)
+
+        C = LazyBandedMatrices.paddeddata(LazyBandedMatrices.invdiagtrav(coefficients(𝐚)))
+        m,n = size(C)
+        using RecurrenceRelationshipArrays
+        X_T = jacobimatrix(T)
+        X_U = jacobimatrix(U)
+        cfs = [Clenshaw(C[1:m-j+1,j], recurrencecoefficients(T)..., X_T) for j=1:n]
+
+    
+        KR = Block.(1:3)
+        
+        @test (KronTrav(Eye(∞),cfs[1]) + KronTrav(2X_U,cfs[2]) + KronTrav((4X_U^2 - I),cfs[3]))[KR,KR] ≈
+                KronTrav(Eye(3),cfs[1][1:3,1:3]) + KronTrav(2X_U[1:3,1:3],cfs[2][1:3,1:3])+ KronTrav((4X_U^2 - I)[1:3,1:3],cfs[3][1:3,1:3]) ≈ A[KR,KR]
+    end
 end