ClenshawKron implemented

Author: dlfivefifty

src/MultivariateOrthogonalPolynomials.jl

@@ -4,7 +4,7 @@ using QuasiArrays: AbstractVector
 using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, BlockArrays, DomainSets,
       QuasiArrays, StaticArrays, ContinuumArrays, InfiniteArrays, InfiniteLinearAlgebra,
       LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
-      HarmonicOrthogonalPolynomials
+      HarmonicOrthogonalPolynomials, RecurrenceRelationships
 
 import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary
 import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle
@@ -21,7 +21,7 @@ import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout, PaddedLayo
 import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes, DiagTravLayout, invdiagtrav, ApplyBandedBlockBandedLayout
 import InfiniteArrays: InfiniteCardinal, OneToInf
 
-import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix
+import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw
 import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
                                           PartialDerivative, AngularMomentum, BlockOneTo, BlockRange1, interlace,
                                           MultivariateOPLayout, MAX_PLOT_BLOCKS
@@ -36,6 +36,7 @@ export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        grammatrix, oneto
 
 include("ModalInterlace.jl")
+include("clenshawkron.jl")
 include("rect.jl")
 include("disk.jl")
 include("rectdisk.jl")

src/clenshawkron.jl

@@ -0,0 +1,51 @@
+
+"""
+    ClenshawKron(coefficientmatrix, (A₁,A₂), (B₁,B₂), (C₁,C₂), (X₁,X₂)))
+
+represents the multiplication operator of a tensor product of two orthogonal polynomials.  That is,
+if `f(x,y) = P(x)*F*Q(y)'` and `a(x,y) = R(x)*A*S(y)'` then
+
+    M = ClenshawKron(A, tuple.(recurrencecoefficients(R), recurrencecoefficients(S)), (jacobimatrix(P), jacobimatrix(Q))
+    M * KronTrav(F)
+
+Gives the coefficients of `a(x,y)*f(x,y)`.
+"""
+struct ClenshawKron{T, Coefs<:AbstractMatrix{T}, Rec<:NTuple{2,NTuple{3,AbstractVector}}, Jac<:NTuple{2,AbstractMatrix}} <: AbstractBandedBlockBandedMatrix{T}
+    c::Coefs
+    recurrencecoeffients::Rec
+    X::Jac
+end
+
+
+copy(M::ClenshawKron) = M
+size(M::ClenshawKron) = (ℵ₀, ℵ₀)
+axes(M::ClenshawKron) = (blockedrange(oneto(∞)),blockedrange(oneto(∞)))
+blockbandwidths(M::ClenshawKron) = (size(M.c,1)-1,size(M.c,1)-1)
+subblockbandwidths(M::ClenshawKron) = (size(M.c,2)-1,size(M.c,2)-1)
+
+function square_getindex(M::ClenshawKron, N::Block{1})
+    # Consider P(x) = L^1_x \ 𝐞_0 
+    # So that if a(x,y) = P(x)*c*Q(y)' then we have
+    # a(X,Y) = 𝐞_0' * inv(L^1_X') * c * Q(Y)'
+    # So we first apply 1D clenshaw to each column
+
+    (A₁,B₁,C₁), (A₂,B₂,C₂) = M.recurrencecoeffients
+    X,Y = M.X
+    m,n = size(M.c)
+    @assert m == n
+    # Apply Clenshaw to each column
+    g = (a,b,N) -> LazyBandedMatrices.krontrav(a[1:N,1:N], b[1:N,1:N])
+    cols = [Clenshaw(M.c[1:m-j+1,j], A₁, B₁, C₁, X) for j=1:n]
+
+    M = m-2+Int(N) # over sample
+    Q_Y = forwardrecurrence(n, A₂, B₂, C₂, Y[1:M,1:M])
+    +(broadcast((a,b,N) -> LazyBandedMatrices.krontrav(a[1:N,1:N], b[1:N,1:N]), Q_Y, cols, Int(N))...)
+end
+
+getindex(M::ClenshawKron, KR::BlockRange{1}, JR::BlockRange{1}) = square_getindex(M, max(maximum(KR), maximum(JR)))[KR,JR]
+getindex(M::ClenshawKron, K::Block{1}, J::Block{1}) = square_getindex(M, max(K, J))[K,J]
+getindex(M::ClenshawKron, Kk::BlockIndex{1}, Jj::BlockIndex{1}) = M[block(Kk), block(Jj)][blockindex(Kk), blockindex(Jj)]
+getindex(M::ClenshawKron, k::Int, j::Int) = M[findblockindex(axes(M,1),k), findblockindex(axes(M,2),j)]
+
+Base.array_summary(io::IO, C::ClenshawKron{T}, inds::Tuple{Vararg{OneToInf{Int}}}) where T =
+    print(io, Base.dims2string(length.(inds)), " ClenshawKron{$T} with $(size(C.c)) polynomial")

test/test_rect.jl

@@ -1,7 +1,8 @@
 using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, LinearAlgebra, BlockArrays, FillArrays, Base64, LazyBandedMatrices, Test
 using ClassicalOrthogonalPolynomials: expand, coefficients, recurrencecoefficients
-using MultivariateOrthogonalPolynomials: weaklaplacian
+using MultivariateOrthogonalPolynomials: weaklaplacian, ClenshawKron
 using ContinuumArrays: plotgridvalues
+using Base: oneto
 
 @testset "RectPolynomial" begin
     @testset "Evaluation" begin
@@ -168,27 +169,19 @@ using ContinuumArrays: plotgridvalues
 
         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, RecurrenceRelationships
-        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]
+        A = ClenshawKron(C, (recurrencecoefficients(T), recurrencecoefficients(U)), (jacobimatrix(T), jacobimatrix(U)))
 
-        g = (a,b,N) -> LazyBandedMatrices.krontrav(a[1:N,1:N], b[1:N,1:N])
-        N = 10
-        M = m-2+N # over sample
-        @time U_X = forwardrecurrence(length(cfs), recurrencecoefficients(U)..., X_U[1:M,1:M]);
-        @time Ks = broadcast(g, U_X, cfs, N);
-        @time A_N = +(Ks...)
-        @test A_N ≈ A[Block.(1:N),Block.(1:N)]
+
+        Ã = a(X,Y)
+        for (k,j) in ((Block.(oneto(5)),Block.(oneto(5))), Block.(oneto(5)),Block.(oneto(6)), (Block(2), Block(3)), (4,5),
+                    (Block(2)[2], Block(3)[3]), (Block(2)[2], Block(3)))
+            @test A[k,j] ≈ Ã[k,j]
+        end
+
+        @test A[Block(1,2)] ≈ Ã[Block(1,2)]
+        @test A[Block(1,2)][1,2] ≈ Ã[Block(1,2)[1,2]]
     end
 end