Faster indexing for Ultraspherical(0.5) ⟷ Legendre (#259)

* Faster indexing for Ultraspherical(0.5)<->Legendre

* short-circuit matching spaces in conversion

Author: jishnub

Project.toml

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

src/Spaces/Jacobi/JacobiOperators.jl

@@ -379,9 +379,14 @@ end
 function getindex(C::ConcreteConversion{<:Ultraspherical,<:Jacobi,T}, k::Integer, j::Integer) where {T}
     if j==k
         S=rangespace(C)
-        jp=jacobip(T,k-1,S.a,S.b,one(T))
-        um=strictconvert(Operator{T}, Evaluation(setcanonicaldomain(domainspace(C)),rightendpoint,0))[k]::T
-        (um/jp)::T
+        U=domainspace(C)
+        if S.a == S.b == 0 && isequalhalf(order(U))
+            oneunit(T)
+        else
+            jp=jacobip(T,k-1,S.a,S.b,one(T))
+            um=strictconvert(Operator{T}, Evaluation(setcanonicaldomain(U),RightEndPoint,0))[k]::T
+            (um/jp)::T
+        end
     else
         zero(T)
     end
@@ -392,12 +397,18 @@ function BandedMatrix(S::SubOperator{T,ConcreteConversion{US,J,T},NTuple{2,UnitR
     kr,jr = parentindices(S)
     k=(kr ∩ jr)
 
-    sp=rangespace(parent(S))
-    jp=jacobip(T,k.-1,sp.a,sp.b,one(T))
-    um=Evaluation(T,setcanonicaldomain(domainspace(parent(S))),rightendpoint,0)[k]
-    vals = um./jp
+    C = parent(S)
+    sp = rangespace(C)
+    U = domainspace(C)
+    if sp.a == 0 && sp.b == 0 && isequalhalf(order(U))
+        ret[band(bandshift(S))] .= oneunit(T)
+    else
+        jp=jacobip(T,k.-1,sp.a,sp.b,one(T))
+        um=Evaluation(T,setcanonicaldomain(U),RightEndPoint,0)[k]
+        vals = um./jp
 
-    ret[band(bandshift(S))] = vals
+        ret[band(bandshift(S))] = vals
+    end
     ret
 end
 
@@ -406,9 +417,14 @@ end
 function getindex(C::ConcreteConversion{<:Jacobi,<:Ultraspherical,T}, k::Integer, j::Integer) where {T}
     if j==k
         S=domainspace(C)
-        jp=jacobip(T,k-1,S.a,S.b,one(T))
-        um=Evaluation(T,setcanonicaldomain(rangespace(C)),rightendpoint,0)[k]
-        jp/um::T
+        U = rangespace(C)
+        if S.a == S.b == 0 && isequalhalf(order(U))
+            oneunit(T)
+        else
+            jp=jacobip(T,k-1,S.a,S.b,one(T))
+            um=Evaluation(T,setcanonicaldomain(rangespace(C)),RightEndPoint,0)[k]
+            jp/um::T
+        end
     else
         zero(T)
     end
@@ -419,12 +435,17 @@ function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,US,T},NTuple{2,UnitR
     kr,jr = parentindices(S)
     k=(kr ∩ jr)
 
-    sp=domainspace(parent(S))
-    jp=jacobip(T,k.-1,sp.a,sp.b,one(T))
-    um=Evaluation(T,setcanonicaldomain(rangespace(parent(S))),rightendpoint,0)[k]
-    vals = jp./um
-
-    ret[band(bandshift(S))] = vals
+    C = parent(S)
+    sp=domainspace(C)
+    U = rangespace(C)
+    if sp.a == 0 && sp.b == 0 && isequalhalf(order(U))
+        ret[band(bandshift(S))] .= oneunit(T)
+    else
+        jp=jacobip(T,k.-1,sp.a,sp.b,one(T))
+        um=Evaluation(T,setcanonicaldomain(U),RightEndPoint,0)[k]
+        vals = jp./um
+        ret[band(bandshift(S))] = vals
+    end
     ret
 end
 

src/Spaces/PolynomialSpace.jl

@@ -466,6 +466,17 @@ function Conversion(L::S, M::NormalizedPolynomialSpace{S}) where S<:PolynomialSp
     end
 end
 
+function Conversion(L::NormalizedPolynomialSpace{<:PolynomialSpace},
+        M::NormalizedPolynomialSpace{<:PolynomialSpace})
+
+    L == M && return Conversion(L)
+
+    C1 = ConcreteConversion(L, canonicalspace(L))
+    C2 = Conversion(canonicalspace(L), canonicalspace(M))
+    C3 = ConcreteConversion(canonicalspace(M), M)
+    ConversionWrapper(C3 * C2 * C1, L, M)
+end
+
 function Fun(::typeof(identity), S::NormalizedPolynomialSpace)
     C = canonicalspace(S)
     f = Fun(identity, C)

test/JacobiTest.jl

@@ -124,7 +124,8 @@ include("testutils.jl")
         end
 
         @testset for (S1, S2) in ((Legendre(), Jacobi(-0.5, -0.5)), (Jacobi(-0.5, -0.5), Legendre()),
-                (Legendre(), Jacobi(1.5, 1.5)))
+                (Legendre(), Jacobi(1.5, 1.5)), (Legendre(), Ultraspherical(0.5)),
+                (Ultraspherical(0.5), Legendre()))
             @test Conversion(S1, S2) * Fun(x->x^4, S1) ≈ Fun(x->x^4, S2)
         end
 
@@ -242,6 +243,39 @@ include("testutils.jl")
             @test domainspace(Y) == Legendre()
             @test Y * Fun(Legendre()) ≈ Fun(x->x^2, Legendre())
         end
+
+        @testset "getindex for ConcreteConversion" begin
+            for m in (0,1)
+                C = Conversion(Jacobi(m,m), Ultraspherical(m+0.5))
+                @test isdiag(C)
+                B = C[1:10, 1:10]
+                for i in 1:10
+                    @test C[i,i] ≈ B[i,i]
+                end
+
+                C = Conversion(Ultraspherical(m+0.5), Jacobi(m,m))
+                @test isdiag(C)
+                B = C[1:10, 1:10]
+                for i in 1:10
+                    @test C[i,i] ≈ B[i,i]
+                end
+            end
+        end
+
+        @testset "NormalizedJacobi and NormalizedUltraspherical" begin
+            for m in (0,1)
+                C = Conversion(NormalizedUltraspherical(m + 0.5), NormalizedJacobi(m,m))
+                @test isdiag(C)
+                @test C[1:10, 1:10] ≈ I
+            end
+        end
+        @testset "Normalized and Unnormalized" begin
+            C = Conversion(NormalizedUltraspherical(0.5), Legendre())
+            @test C * Fun(x->x^4, NormalizedUltraspherical(0.5)) ≈ Fun(x->x^4, Legendre())
+
+            C = Conversion(Legendre(), NormalizedChebyshev())
+            @test C * Fun(x->x^4, Legendre()) ≈ Fun(x->x^4, NormalizedChebyshev())
+        end
     end
 
     @testset "inplace transform" begin