Specialize isdiag for certain ConversionWrappers (#287)

* Specialize isdiag for certain ConversionWrappers

* isdiag for Jacobi ConversionWrapper

* Bump version to v0.6.42

* Conversions between Jacobi and Ultraspherical

* Specialize israggedbelow, fix Ultraspheriacl ConcreteConversion indexing

* Specialize ConcreteConversion indexing for Union{Integer, HalfOddInteger}

* Fix ambiguity

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.41"
+version = "0.6.42"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"

src/ApproxFunOrthogonalPolynomials.jl

@@ -56,7 +56,8 @@ import ApproxFunBase: Fun, SubSpace, WeightSpace, NoSpace, HeavisideSpace,
                     isleftendpoint, isrightendpoint, evaluation_point, boundaryfn, ldiffbc, rdiffbc,
                     LeftEndPoint, RightEndPoint, normalizedspace, promotedomainspace,
                     bandmatrices_eigen, SymmetricEigensystem, SkewSymmetricEigensystem,
-                    mean # differs from Statistics.mean after https://github.com/JuliaApproximation/ApproxFunBase.jl/pull/506
+                    mean, # differs from Statistics.mean after https://github.com/JuliaApproximation/ApproxFunBase.jl/pull/506
+                    israggedbelow
 
 import DomainSets: Domain, indomain, UnionDomain, FullSpace, Point,
             Interval, ChebyshevInterval, boundary, rightendpoint, leftendpoint,

src/Spaces/Jacobi/JacobiOperators.jl

@@ -192,7 +192,14 @@ end
 
 bandwidths(::ConcreteConversion{<:Jacobi,<:Jacobi}) = (0,1)
 
-
+# These methods help with constant propagating the bandwidths of a KroneckerOperator
+isdiag(::ConversionWrapper{<:Chebyshev,<:Jacobi}) = false
+isdiag(::ConversionWrapper{<:Jacobi,<:Jacobi,<:Any,<:ConstantOperator}) = true
+isdiag(::ConversionWrapper{<:Union{Ultraspherical,Jacobi},<:Union{Ultraspherical,Jacobi}}) = false
+
+# These operators are banded and UpperTriangular
+israggedbelow(::ConversionWrapper{<:Chebyshev,<:Jacobi}) = true
+israggedbelow(::ConversionWrapper{<:Union{Ultraspherical,Jacobi},<:Union{Ultraspherical,Jacobi}}) = true
 
 function getindex(C::ConcreteConversion{<:Jacobi,<:Jacobi,T},k::Integer,j::Integer) where {T}
     L=C.domainspace

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -213,10 +213,10 @@ function getindex(M::ConcreteConversion{<:Chebyshev,U,T},
 end
 
 
-function getindex(M::ConcreteConversion{U1,U2,T},
-        k::Integer,j::Integer) where {DD,RR,
-            U1<:Ultraspherical{<:Integer,DD,RR},
-            U2<:Ultraspherical{<:Integer,DD,RR},T}
+function _getindex(M::ConcreteConversion{U1,U2,T}, k, j) where {DD,RR,
+            OT<:Union{Integer, HalfOddInteger},
+            U1<:Ultraspherical{<:OT,DD,RR},
+            U2<:Ultraspherical{<:OT,DD,RR},T}
     #  we can assume that λ==m+1
     λ=order(rangespace(M))
     c=λ-one(T)  # this supports big types
@@ -229,6 +229,17 @@ function getindex(M::ConcreteConversion{U1,U2,T},
     end
 end
 
+function getindex(M::ConcreteConversion{U1,U2,T}, k::Integer, j::Integer) where {DD,RR,
+            U1<:Ultraspherical{<:Integer,DD,RR},
+            U2<:Ultraspherical{<:Integer,DD,RR},T}
+    _getindex(M, k, j)
+end
+function getindex(M::ConcreteConversion{U1,U2,T}, k::Integer, j::Integer) where {DD,RR,
+            U1<:Ultraspherical{<:HalfOddInteger,DD,RR},
+            U2<:Ultraspherical{<:HalfOddInteger,DD,RR},T}
+    _getindex(M, k, j)
+end
+
 function getindex(M::ConcreteConversion{U1,U2,T},
         k::Integer,j::Integer) where {DD,RR,
             U1<:Ultraspherical{<:Any,DD,RR},
@@ -273,6 +284,13 @@ isdiag(::ConcreteConversion{<:Chebyshev,<:Ultraspherical}) = false
 isdiag(::ConcreteConversion{<:Ultraspherical,<:Chebyshev}) = false
 isdiag(::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) = false
 
+# These methods help with constant propagating the bandwidths of a KroneckerOperator
+isdiag(::ConversionWrapper{<:Chebyshev,<:Ultraspherical}) = false
+isdiag(::ConversionWrapper{<:Ultraspherical,<:Ultraspherical, <:Any, <:ConstantOperator}) = true
+
+# These operators are banded and UpperTriangular
+israggedbelow(::ConversionWrapper{<:Chebyshev,<:Ultraspherical}) = true
+
 ## coefficients
 
 # return the space that has banded Conversion to the other
@@ -381,7 +399,7 @@ function getindex(M::ConcreteConversion{Ultraspherical{LT,DD,RR},
             zero(T)
         end
     else
-        error("Not implemented")
+        error("Not implemented for $M")
     end
 end
 

test/JacobiTest.jl

@@ -791,6 +791,63 @@ include("testutils.jl")
         @test ApproxFunBase.hasconversion(NormalizedChebyshev()*Legendre(), Chebyshev()*Legendre())
         @test ApproxFunBase.hasconversion(NormalizedChebyshev()*NormalizedLegendre(), Chebyshev()*Legendre())
         @test ApproxFunBase.hasconversion(Chebyshev()*Legendre(), NormalizedChebyshev()*NormalizedLegendre())
+
+        @testset "isbanded and israggedbelow for KroneckerOperator" begin
+            function testspaces(d1, d2, r1, r2, res = Val(false))
+                K = (Operator(I, d1) ⊗ Operator(I, d2)) → (r1 * r2)
+                @test (@inferred (K -> Val(isbanded(K)))(K)) == res
+                @test (@inferred (K -> Val(ApproxFunBase.israggedbelow(K)))(K)) == Val(true)
+                f = Fun((x,y)->x*y, d1*d2)
+                g = Fun((x,y)->x*y, r1*r2)
+                @test K * f ≈ g
+            end
+
+            d1, r1 = Chebyshev(), Legendre()
+            d2, r2 = Chebyshev(), Chebyshev()
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Legendre(), Chebyshev()
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Chebyshev(), Jacobi(1,1)
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Chebyshev(), Jacobi(Ultraspherical(1))
+            d2, r2 = Chebyshev(), Chebyshev()
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Chebyshev(), Jacobi(Ultraspherical(2))
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Jacobi(Ultraspherical(1)), Jacobi(Ultraspherical(1))
+            testspaces(d1, d2, r1, r2, Val(true))
+
+            d1, r1 = Jacobi(Ultraspherical(1)), Jacobi(Ultraspherical(2))
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Jacobi(Ultraspherical(1)), Ultraspherical(2)
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Ultraspherical(1), Jacobi(Ultraspherical(2))
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Jacobi(Ultraspherical(1)), Jacobi(Ultraspherical(3))
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Jacobi(Ultraspherical(1)), Jacobi(2,2)
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Ultraspherical(1), Jacobi(2,2)
+            testspaces(d1, d2, r1, r2)
+        end
+
+        d1, r1 = Legendre(), Jacobi(2,2)
+        d2, r2 = Chebyshev(), Chebyshev()
+        K = (Operator(I, d1) ⊗ Operator(I, d2)) → (r1 ⊗ r2)
+        M = K[1:20, 1:20]
+        for ind in CartesianIndices(M)
+            @test K[ind] ≈ M[ind]
+        end
     end
 end
 

test/UltrasphericalTest.jl

@@ -2,12 +2,13 @@ module UltrasphericalTest
 
 using ApproxFunOrthogonalPolynomials
 using ApproxFunBase
-using Test
-using SpecialFunctions
+using BandedMatrices
+using HalfIntegers
 using LinearAlgebra
-using Static
 using OddEvenIntegers
-using HalfIntegers
+using SpecialFunctions
+using Static
+using Test
 
 include("testutils.jl")
 
@@ -187,6 +188,45 @@ include("testutils.jl")
                 end
             end
         end
+
+        @testset "isbanded for KroneckerOperator" begin
+            function testspaces(d1, d2, r1, r2, res = Val(false))
+                K = (Operator(I, d1) ⊗ Operator(I, d2)) → (r1 * r2)
+                @test (@inferred (K -> Val(isbanded(K)))(K)) == res
+                @test (@inferred (K -> Val(ApproxFunBase.israggedbelow(K)))(K)) == Val(true)
+                f = Fun((x,y)->x*y, d1*d2)
+                g = Fun((x,y)->x*y, r1*r2)
+                @test K * f ≈ g
+            end
+            d1, r1 = Chebyshev(), Ultraspherical(1)
+            d2, r2 = Chebyshev(), Chebyshev()
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Chebyshev(), Ultraspherical(2)
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Ultraspherical(1), Ultraspherical(1)
+            testspaces(d1, d2, r1, r2, Val(true))
+
+            d1, r1 = Ultraspherical(1), Ultraspherical(2)
+            testspaces(d1, d2, r1, r2)
+
+            d1, r1 = Ultraspherical(1), Ultraspherical(3)
+            testspaces(d1, d2, r1, r2)
+        end
+
+        @testset "ConcreteConversion indexing" begin
+            for (s1, s2) in ((Ultraspherical(Legendre()), Ultraspherical(Jacobi(1,1))),
+                            (Ultraspherical(1), Ultraspherical(2)),
+                            (Ultraspherical(static(1)), Ultraspherical(static(2))),
+                            )
+                C = Conversion(s1, s2)
+                M = C[1:4, 1:4]
+                for ind in CartesianIndices(M)
+                    @test C[ind] ≈ M[ind]
+                end
+            end
+        end
     end
 
     @testset "Normalized space" begin