Concrete derivatives for Ultraspherical and Jacobi

Author: jishnub

src/Spaces/Chebyshev/Chebyshev.jl

@@ -26,7 +26,7 @@ const NormalizedChebyshev{D<:Domain,R} = NormalizedPolynomialSpace{Chebyshev{D,
 NormalizedChebyshev() = NormalizedPolynomialSpace(Chebyshev())
 NormalizedChebyshev(d) = NormalizedPolynomialSpace(Chebyshev(d))
 
-function Base.getproperty(S::Chebyshev{<:Any,<:Any},v::Symbol)
+function getproperty(S::Chebyshev, v::Symbol)
     if v==:b || v==:a
         -0.5
     else

src/Spaces/Jacobi/JacobiOperators.jl

@@ -1,7 +1,7 @@
 ## Derivative
 
 # specialize Derivative so that this is type-inferred even without constant propagation
-Derivative(J::Jacobi) = ConcreteDerivative(J,1)
+Derivative(J::MaybeNormalized{<:Jacobi}) = ConcreteDerivative(J,1)
 @inline function _Derivative(J::Jacobi, k::Number)
     assert_integer(k)
     if k==1
@@ -20,9 +20,9 @@ else
 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
+rangespace(D::ConcreteDerivative{<:MaybeNormalized{<:Jacobi}}) = Jacobi(D.space.b+D.order,D.space.a+D.order,domain(D))
+bandwidths(D::ConcreteDerivative{<:MaybeNormalized{<:Jacobi}}) = -D.order,D.order
+isdiag(D::ConcreteDerivative{<:MaybeNormalized{<: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))

src/Spaces/PolynomialSpace.jl

@@ -430,6 +430,11 @@ end
 
 @containsconstants NormalizedPolynomialSpace
 
+@inline function getproperty(N::NormalizedPolynomialSpace, v::Symbol)
+    ((v == :a) || (v == :b)) && return getproperty(getfield(N, :space), v)
+    getfield(N, v)
+end
+
 domain(S::NormalizedPolynomialSpace) = domain(S.space)
 canonicalspace(S::NormalizedPolynomialSpace) = S.space
 setdomain(NS::NormalizedPolynomialSpace, d::Domain) = NormalizedPolynomialSpace(setdomain(canonicalspace(NS), d))
@@ -505,8 +510,8 @@ end
 
 # Conversion
 
-bandwidths(C::ConcreteConversion{NormalizedPolynomialSpace{S,D,R},S}) where {S,D,R} = (0, 0)
-bandwidths(C::ConcreteConversion{S,NormalizedPolynomialSpace{S,D,R}}) where {S,D,R} = (0, 0)
+bandwidths(::ConcreteConversion{<:NormalizedPolynomialSpace{S},S}) where {S} = (0, 0)
+bandwidths(::ConcreteConversion{S,<:NormalizedPolynomialSpace{S}}) where {S} = (0, 0)
 
 function getindex(C::ConcreteConversion{NormalizedPolynomialSpace{S,D,R},S,T},k::Integer,j::Integer) where {S,D<:IntervalOrSegment,R,T}
     if j==k

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -57,7 +57,7 @@ end
 #Derivative(d::IntervalOrSegment)=Derivative(1,d)
 
 
-function Derivative(sp::Ultraspherical{LT,DD}, m::Number) where {LT,DD<:IntervalOrSegment}
+function Derivative(sp::MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}, m::Number)
     assert_integer(m)
     ConcreteDerivative(sp,m)
 end
@@ -75,15 +75,15 @@ function Integral(sp::Ultraspherical{<:Any,<:IntervalOrSegment}, m::Number)
 end
 
 
-rangespace(D::ConcreteDerivative{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} =
+rangespace(D::ConcreteDerivative{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) =
     Ultraspherical(order(domainspace(D))+D.order,domain(D))
 
-bandwidths(D::ConcreteDerivative{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} = -D.order,D.order
-bandwidths(D::ConcreteIntegral{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} = D.order,-D.order
-Base.stride(D::ConcreteDerivative{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} = D.order
+bandwidths(D::ConcreteDerivative{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = -D.order,D.order
+bandwidths(D::ConcreteIntegral{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = D.order,-D.order
+Base.stride(D::ConcreteDerivative{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = D.order
 
-isdiag(D::ConcreteDerivative{<:Ultraspherical{<:Any,<:IntervalOrSegment}}) = false
-isdiag(D::ConcreteIntegral{<:Ultraspherical{<:Any,<:IntervalOrSegment}}) = false
+isdiag(D::ConcreteDerivative{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = false
+isdiag(D::ConcreteIntegral{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = false
 
 function getindex(D::ConcreteDerivative{<:Ultraspherical{TT,DD},K,T},
                k::Integer,j::Integer) where {TT,DD<:IntervalOrSegment,K,T}

src/fastops.jl

@@ -72,8 +72,8 @@ end
 
 
 
-function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Chebyshev{DD,RR},K,T},
-                                                     NTuple{2,UnitRange{Int}}}) where {T,K,DD,RR}
+function BandedMatrix(S::SubOperator{T,<:ConcreteDerivative{<:Chebyshev,<:Any,T},
+                                                     NTuple{2,UnitRange{Int}}}) where {T}
 
     n,m = size(S)
     ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
@@ -98,8 +98,8 @@ function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Chebyshev{DD,RR},K,T},
 end
 
 
-function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Ultraspherical{LT,DD,RR},K,T},
-                                                  NTuple{2,UnitRange{Int}}}) where {T,K,DD,RR,LT}
+function BandedMatrix(S::SubOperator{T,<:ConcreteDerivative{<:Ultraspherical,T},
+                                                  NTuple{2,UnitRange{Int}}}) where {T}
     n,m = size(S)
     ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
     kr,jr = parentindices(S)
@@ -116,3 +116,18 @@ function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Ultraspherical{LT,DD,R
     ret
 end
 
+function BandedMatrix(S::SubOperator{T,<:ConcreteDerivative{<:NormalizedPolynomialSpace,<:Any,T},
+                                                     NTuple{2,UnitRange{Int}}}) where {T}
+
+    D = parent(S)
+    sp = domainspace(D)
+    csp = canonicalspace(sp)
+    ind1, ind2 = parentindices(S)
+    B = BandedMatrix(view(ConcreteDerivative(csp, D.order), ind1, ind2))
+    data = BandedMatrices.bandeddata(B)
+    C = ConcreteConversion(sp, csp)
+    for (ind, col) in enumerate(ind2)
+        @views data[:, ind] .*= C[col, col]
+    end
+    B
+end