isdiag for known derivatives (#183)

* isdiag for known derivatives

* isdiag for NormalizedPolynomialSpace

* Concreteint Jacobi and Concreteconv Ultrasp

* Version bump to v0.6.7

* Integral tests

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.6"
+version = "0.6.7"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
@@ -20,7 +20,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ApproxFunBase = "0.7.65"
+ApproxFunBase = "0.7.67"
 ApproxFunBaseTest = "0.1"
 Aqua = "0.5"
 BandedMatrices = "0.16, 0.17"

src/ApproxFunOrthogonalPolynomials.jl

@@ -79,6 +79,8 @@ import SpecialFunctions: erfcx, dawson,
 
 using StaticArrays: SVector
 
+import LinearAlgebra: isdiag
+
 points(d::IntervalOrSegmentDomain{T},n::Integer) where {T} =
     fromcanonical.(Ref(d), chebyshevpoints(float(real(eltype(T))), n))  # eltype to handle point
 bary(v::AbstractVector{Float64},d::IntervalOrSegmentDomain,x::Float64) = bary(v,tocanonical(d,x))

src/Spaces/Chebyshev/ChebyshevOperators.jl

@@ -226,6 +226,8 @@ rangespace(D::ConcreteDerivative{Chebyshev{DD,RR}}) where {DD<:IntervalOrSegment
 bandwidths(D::ConcreteDerivative{Chebyshev{DD,RR}}) where {DD<:IntervalOrSegment,RR} = -D.order,D.order
 Base.stride(D::ConcreteDerivative{Chebyshev{DD,RR}}) where {DD<:IntervalOrSegment,RR} = D.order
 
+isdiag(D::ConcreteDerivative{<:Chebyshev{<:IntervalOrSegment}}) = false
+
 function getindex(D::ConcreteDerivative{Chebyshev{DD,RR},K,T},k::Integer,j::Integer) where {DD<:IntervalOrSegment,RR,K,T}
     m=D.order
     d=domain(D)

src/Spaces/Hermite/Hermite.jl

@@ -45,6 +45,8 @@ rangespace(D::ConcreteDerivative{H}) where {H<:Hermite} = domainspace(D)
 getindex(D::ConcreteDerivative{H},k::Integer,j::Integer) where {H<:Hermite} =
         j==k+D.order ? one(eltype(D))*2^D.order*pochhammer(k,D.order) : zero(eltype(D))
 
+isdiag(D::ConcreteDerivative{<:Hermite}) = false
+
 function hermitep(r::AbstractRange,x::Number)
     n = r[end] + 1
 

src/Spaces/Jacobi/JacobiOperators.jl

@@ -20,6 +20,7 @@ end
 
 rangespace(D::ConcreteDerivative{<:Jacobi}) = Jacobi(D.space.b+D.order,D.space.a+D.order,domain(D))
 bandwidths(D::ConcreteDerivative{<:Jacobi}) = -D.order,D.order
+isdiag(D::ConcreteDerivative{<:Jacobi}) = false
 
 getindex(T::ConcreteDerivative{<:Jacobi}, k::Integer, j::Integer) =
     j==k+1 ? eltype(T)((k+1+T.space.a+T.space.b)/complexlength(domain(T))) : zero(eltype(T))
@@ -43,6 +44,7 @@ Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::Number,o::Integer) = Concre
 
 function Integral(J::Jacobi,k::Number)
     assert_integer(k)
+    @assert k > 0 "order of integral must be > 0"
     if k > 1
         Q=Integral(J,1)
         IntegralWrapper(TimesOperator(Integral(rangespace(Q),k-1),Q),J,k)
@@ -58,7 +60,8 @@ end
 
 
 rangespace(D::ConcreteIntegral{<:Jacobi}) = Jacobi(D.space.b-D.order,D.space.a-D.order,domain(D))
-bandwidths(D::ConcreteIntegral{<:Jacobi}) = D.order,0
+bandwidths(D::ConcreteIntegral{<:Jacobi}) = D.order,-D.order
+isdiag(D::ConcreteIntegral{<:Jacobi}) = false
 
 function getindex(T::ConcreteIntegral{<:Jacobi}, k::Integer, j::Integer)
     @assert T.order==1

src/Spaces/PolynomialSpace.jl

@@ -456,6 +456,12 @@ hasconversion(a::PolynomialSpace,b::NormalizedPolynomialSpace) = hasconversion(a
 hasconversion(a::NormalizedPolynomialSpace,b::PolynomialSpace) = hasconversion(a.space,b)
 hasconversion(a::NormalizedPolynomialSpace,b::NormalizedPolynomialSpace) = hasconversion(a.space,b)
 
+function isdiag(D::DerivativeWrapper{<:Operator, <:NormalizedPolynomialSpace})
+    sp = D.space
+    csp = _stripnorm(sp)
+    isdiag(Derivative(csp, D.order))
+end
+
 # Tensor products of normalized and unnormalized spaces may have banded conversions defined
 # A banded conversion exists in special cases, where both conversion operators are diagonal
 _stripnorm(N::NormalizedPolynomialSpace) = canonicalspace(N)

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -77,6 +77,8 @@ bandwidths(D::ConcreteDerivative{Ultraspherical{LT,DD,RR}}) where {LT,DD<:Interv
 bandwidths(D::ConcreteIntegral{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} = D.order,-D.order
 Base.stride(D::ConcreteDerivative{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} = D.order
 
+isdiag(D::ConcreteDerivative{<:Ultraspherical{<:Any,<:IntervalOrSegment}}) = false
+isdiag(D::ConcreteIntegral{<:Ultraspherical{<:Any,<:IntervalOrSegment}}) = false
 
 function getindex(D::ConcreteDerivative{Ultraspherical{TT,DD,RR},K,T},
                k::Integer,j::Integer) where {TT,DD<:IntervalOrSegment,RR,K,T}
@@ -230,6 +232,9 @@ bandwidths(C::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) =
 Base.stride(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical{Int}}) = 2
 Base.stride(C::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) = 2
 
+isdiag(::ConcreteConversion{<:Chebyshev,<:Ultraspherical}) = false
+isdiag(::ConcreteConversion{<:Ultraspherical,<:Chebyshev}) = false
+isdiag(::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) = false
 
 ## coefficients
 

test/HermiteTest.jl

@@ -68,6 +68,7 @@ using ApproxFunBaseTest: testbandedoperator
     @testset "Derivative" begin
         D = Derivative(Hermite())
         testbandedoperator(D)
+        @test !isdiag(D)
 
         f = Fun( x-> x + x^2, Hermite())
         g = D * f

test/JacobiTest.jl

@@ -180,6 +180,9 @@ using Static
         # only one band should be populated
         @test bandwidths(D, 1) == -bandwidths(D, 2)
 
+        @test !isdiag(Derivative(Legendre()))
+        @test !isdiag(Derivative(NormalizedLegendre()))
+
         @testset for d in [-1..1, 0..1]
             f = Fun(x->x^2, Chebyshev(d))
             C = space(f)
@@ -535,6 +538,20 @@ using Static
         end
     end
 
+    @testset "Integral" begin
+        @testset for sp in (Legendre(), Jacobi(1,1))
+            Ij = Integral(sp, 1)
+            @test !isdiag(Ij)
+            f = Fun(sp)
+            g = Ij * f
+            g = Fun(g, sp)
+            g = g - coefficients(g)[1]
+            gexp = Fun(x->x^2/2, sp)
+            gexp = gexp - coefficients(gexp)[1]
+            @test g ≈ gexp
+        end
+    end
+
     @testset "type inference" begin
         x = @inferred ApproxFunBase.maxspace(Jacobi(1,1), Jacobi(2,2))
         @test x == Jacobi(2,2)

test/UltrasphericalTest.jl

@@ -179,5 +179,17 @@ using LinearAlgebra
         A = @inferred Integral(d)
         x = Derivative() * A * Fun(d)
         @test x(0.2) ≈ 0.2
+
+        @testset for sp in (Ultraspherical(1), Ultraspherical(2), Ultraspherical(0.5))
+            Ij = Integral(sp, 1)
+            @test !isdiag(Ij)
+            f = Fun(sp)
+            g = Ij * f
+            g = Fun(g, sp)
+            g = g - coefficients(g)[1]
+            gexp = Fun(x->x^2/2, sp)
+            gexp = gexp - coefficients(gexp)[1]
+            @test g ≈ gexp
+        end
     end
 end