Skip to content

static integer orders in Ultraspherical spaces #187

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Feb 17, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions src/Spaces/Jacobi/Jacobi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,13 @@ struct Jacobi{D<:Domain,R,T} <: PolynomialSpace{D,R}
domain::D
Jacobi{D,R}(b::T,a::T,d::D) where {D,R,T} = new{D,R,T}(b,a,d)
end
Jacobi(b::T,a::T,d::Domain) where {T} =
Jacobi(b::T,a::T,d::Domain) where {T<:Number} =
Jacobi{typeof(d),promote_type(T,real(prectype(d)))}(b, a, d)
Legendre(domain) = Jacobi(0,0,domain)
Legendre() = Legendre(ChebyshevInterval())
Jacobi(b,a,d::Domain) = Jacobi(promote(dynamic(b), dynamic(a))...,d)
Jacobi(b,a,d) = Jacobi(b,a,Domain(d))
Jacobi(b,a) = Jacobi(b,a,ChebyshevInterval())
Jacobi(b::Number,a::Number,d::Domain) = Jacobi(promote(dynamic(b), dynamic(a))...,d)
Jacobi(b::Number,a::Number,d) = Jacobi(b,a,Domain(d))
Jacobi(b::Number,a::Number) = Jacobi(b,a,ChebyshevInterval())
Jacobi(A::Ultraspherical) = Jacobi(order(A)-0.5,order(A)-0.5,domain(A))
Jacobi(A::Chebyshev) = Jacobi(-0.5,-0.5,domain(A))

Expand Down
4 changes: 2 additions & 2 deletions src/Spaces/Jacobi/JacobiOperators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -153,12 +153,12 @@ function Conversion(L::Jacobi,M::Jacobi)
ConversionWrapper(
TimesOperator(
ConcreteConversion(Jacobi(M.b-1,M.a,dm),M),
Conversion(L,Jacobi(M.b-1,M.a,dm))))
Conversion(L,Jacobi(M.b-static(1),M.a,dm))))
else #if M.a >= L.a+1
ConversionWrapper(
TimesOperator(
ConcreteConversion(Jacobi(M.b,M.a-1,dm),M),
Conversion(L,Jacobi(M.b,M.a-1,dm))))
Conversion(L,Jacobi(M.b,M.a-static(1),dm))))
end
elseif L.a ≈ L.b ≈ 0 && M.a ≈ M.b ≈ 0.5
Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
Expand Down
32 changes: 12 additions & 20 deletions src/Spaces/Ultraspherical/Ultraspherical.jl
Original file line number Diff line number Diff line change
Expand Up @@ -44,47 +44,39 @@ pointscompatible(A::Chebyshev, B::Ultraspherical) = domain(A) == domain(B)

struct UltrasphericalPlan{CT,FT,IP}
chebplan::CT
cheb2legplan::FT
cheb2ultraplan::FT

UltrasphericalPlan{CT,FT}(cp,c2lp,::Val{IP}) where {CT,FT,IP} = new{CT,FT,IP}(cp,c2lp)
end

struct UltrasphericalIPlan{CT,FT,IP}
chebiplan::CT
leg2chebplan::FT
ultra2chebplan::FT

UltrasphericalIPlan{CT,FT}(cp,c2lp,::Val{IP}) where {CT,FT,IP} = new{CT,FT,IP}(cp,c2lp)
end

function UltrasphericalPlan(λ::Number,vals,inplace = Val(false))
if λ == 0.5
cp = ApproxFunBase._plan_transform!!(inplace)(Chebyshev(),vals)
c2lp = plan_cheb2leg(eltype(vals),length(vals))
UltrasphericalPlan{typeof(cp),typeof(c2lp)}(cp,c2lp,inplace)
else
error("Not implemented")
end
cp = ApproxFunBase._plan_transform!!(inplace)(Chebyshev(),vals)
c2lp = plan_cheb2ultra(vals, λ)
UltrasphericalPlan{typeof(cp),typeof(c2lp)}(cp,c2lp,inplace)
end

function UltrasphericalIPlan(λ::Number,cfs,inplace = Val(false))
if λ == 0.5
cp = ApproxFunBase._plan_itransform!!(inplace)(Chebyshev(),cfs)
c2lp=plan_leg2cheb(eltype(cfs),length(cfs))
UltrasphericalIPlan{typeof(cp),typeof(c2lp)}(cp,c2lp,inplace)
else
error("Not implemented")
end
cp = ApproxFunBase._plan_itransform!!(inplace)(Chebyshev(),cfs)
l2cp = plan_ultra2cheb(cfs, λ)
UltrasphericalIPlan{typeof(cp),typeof(l2cp)}(cp,l2cp,inplace)
end

*(UP::UltrasphericalPlan{<:Any,<:Any,false},v::AbstractVector) =
UP.cheb2legplan*(UP.chebplan*v)
UP.cheb2ultraplan*(UP.chebplan*v)
*(UP::UltrasphericalIPlan{<:Any,<:Any,false},v::AbstractVector) =
UP.chebiplan*(UP.leg2chebplan*v)
UP.chebiplan*(UP.ultra2chebplan*v)

*(UP::UltrasphericalPlan{<:Any,<:Any,true},v::AbstractVector) =
lmul!(UP.cheb2legplan, UP.chebplan*v)
lmul!(UP.cheb2ultraplan, UP.chebplan*v)
*(UP::UltrasphericalIPlan{<:Any,<:Any,true},v::AbstractVector) =
UP.chebiplan * lmul!(UP.leg2chebplan, v)
UP.chebiplan * lmul!(UP.ultra2chebplan, v)

plan_transform(sp::Ultraspherical{Int},vals::AbstractVector) = CanonicalTransformPlan(sp,vals)
plan_transform!(sp::Ultraspherical{Int},vals::AbstractVector) = CanonicalTransformPlan(sp,vals,Val(true))
Expand Down
69 changes: 38 additions & 31 deletions src/Spaces/Ultraspherical/UltrasphericalOperators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,14 +19,14 @@ normalization(::Type{T}, sp::Ultraspherical, k::Int) where T = (λ = order(sp);
# these are special cases


Base.stride(M::ConcreteMultiplication{U,V}) where {U<:Chebyshev,V<:Ultraspherical} =
Base.stride(M::ConcreteMultiplication{<:Chebyshev,<:Ultraspherical}) =
stride(M.f)
Base.stride(M::ConcreteMultiplication{U,V}) where {U<:Ultraspherical,V<:Chebyshev} =
Base.stride(M::ConcreteMultiplication{<:Ultraspherical,<:Chebyshev}) =
stride(M.f)
Base.stride(M::ConcreteMultiplication{U,V}) where {U<:Ultraspherical,V<:Ultraspherical} =
Base.stride(M::ConcreteMultiplication{<:Ultraspherical,<:Ultraspherical}) =
stride(M.f)

@inline function _Multiplication(f::Fun{<:Chebyshev}, sp::Ultraspherical{Int})
@inline function _Multiplication(f::Fun{<:Chebyshev}, sp::Ultraspherical{<:Union{Integer,StaticInt}})
if order(sp) == 1
cfs = f.coefficients
MultiplicationWrapper(f,
Expand All @@ -40,10 +40,10 @@ Base.stride(M::ConcreteMultiplication{U,V}) where {U<:Ultraspherical,V<:Ultrasph
end
end
@static if VERSION >= v"1.8"
Base.@constprop aggressive Multiplication(f::Fun{<:Chebyshev}, sp::Ultraspherical{Int}) =
Base.@constprop aggressive Multiplication(f::Fun{<:Chebyshev}, sp::Ultraspherical{<:Union{Integer,StaticInt}}) =
_Multiplication(f, sp)
else
Multiplication(f::Fun{<:Chebyshev}, sp::Ultraspherical{Int}) = _Multiplication(f, sp)
Multiplication(f::Fun{<:Chebyshev}, sp::Ultraspherical{<:Union{Integer,StaticInt}}) = _Multiplication(f, sp)
end


Expand All @@ -70,18 +70,18 @@ function Integral(sp::Ultraspherical{<:Any,<:IntervalOrSegment}, m::Number)
end


rangespace(D::ConcreteDerivative{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} =
rangespace(D::ConcreteDerivative{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} =
Ultraspherical(order(domainspace(D))+D.order,domain(D))

bandwidths(D::ConcreteDerivative{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} = -D.order,D.order
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
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

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}
function getindex(D::ConcreteDerivative{<:Ultraspherical{TT,DD},K,T},
k::Integer,j::Integer) where {TT,DD<:IntervalOrSegment,K,T}
m=D.order
d=domain(D)
λ=order(domainspace(D))
Expand All @@ -96,14 +96,14 @@ end

## Integral

linesum(f::Fun{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} =
linesum(f::Fun{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} =
sum(setcanonicaldomain(f))*arclength(d)/2


rangespace(D::ConcreteIntegral{Ultraspherical{LT,DD,RR}}) where {LT,DD<:IntervalOrSegment,RR} =
rangespace(D::ConcreteIntegral{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} =
order(domainspace(D)) == 1 ? Chebyshev(domain(D)) : Ultraspherical(order(domainspace(D))-D.order,domain(D))

function getindex(Q::ConcreteIntegral{Ultraspherical{LT,DD,RR}},k::Integer,j::Integer) where {LT,DD<:IntervalOrSegment,RR}
function getindex(Q::ConcreteIntegral{<:Ultraspherical{LT,DD}},k::Integer,j::Integer) where {LT,DD<:IntervalOrSegment}
T=eltype(Q)
m=Q.order
d=domain(Q)
Expand Down Expand Up @@ -160,7 +160,7 @@ function Conversion(A::Ultraspherical,B::Ultraspherical)
ConcreteConversion(A,B)
elseif b-a > 1
d=domain(A)
US=Ultraspherical(b-1,d)
US=Ultraspherical(b-static(1),d)
ConversionWrapper(TimesOperator(
ConcreteConversion(US,B), Conversion(A,US)))
else
Expand All @@ -171,8 +171,8 @@ end
maxspace_rule(A::Ultraspherical,B::Ultraspherical) = order(A) > order(B) ? A : B


function getindex(M::ConcreteConversion{C,Ultraspherical{Int,DD,RR},T},
k::Integer,j::Integer) where {DD,RR,C<:Chebyshev,T}
function getindex(M::ConcreteConversion{<:Chebyshev,U,T},
k::Integer,j::Integer) where {T, U<:Ultraspherical{<:Union{Integer, StaticInt}}}
# order must be 1
if k==j==1
one(T)
Expand All @@ -186,8 +186,10 @@ function getindex(M::ConcreteConversion{C,Ultraspherical{Int,DD,RR},T},
end


function getindex(M::ConcreteConversion{Ultraspherical{Int,DD,RR},Ultraspherical{Int,DD,RR},T},
k::Integer,j::Integer) where {DD,RR,T}
function getindex(M::ConcreteConversion{U1,U2,T},
k::Integer,j::Integer) where {DD,RR,
U1<:Ultraspherical{<:Union{Integer, StaticInt},DD,RR},
U2<:Ultraspherical{<:Union{Integer, StaticInt},DD,RR},T}
# we can assume that λ==m+1
λ=order(rangespace(M))
c=λ-one(T) # this supports big types
Expand All @@ -200,8 +202,10 @@ function getindex(M::ConcreteConversion{Ultraspherical{Int,DD,RR},Ultraspherical
end
end

function getindex(M::ConcreteConversion{Ultraspherical{LT,DD,RR},Ultraspherical{LT,DD,RR},T},
k::Integer,j::Integer) where {LT,DD,RR,T}
function getindex(M::ConcreteConversion{U1,U2,T},
k::Integer,j::Integer) where {DD,RR,
U1<:Ultraspherical{<:Any,DD,RR},
U2<:Ultraspherical{<:Any,DD,RR},T}
λ=order(rangespace(M))
if order(domainspace(M))+1==λ
c=λ-one(T) # this supports big types
Expand All @@ -218,8 +222,8 @@ function getindex(M::ConcreteConversion{Ultraspherical{LT,DD,RR},Ultraspherical{
end


bandwidths(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical{Int}}) = 0,2 # order == 1
bandwidths(C::ConcreteConversion{<:Ultraspherical{Int},<:Ultraspherical{Int}}) = 0,2
bandwidths(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical{<:Union{Integer,StaticInt}}}) = 0,2 # order == 1
bandwidths(C::ConcreteConversion{<:Ultraspherical{<:Union{Integer,StaticInt}},<:Ultraspherical{<:Union{Integer,StaticInt}}}) = 0,2

bandwidths(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical}) =
0,order(rangespace(C))==1 ? 2 : ℵ₀
Expand All @@ -229,7 +233,7 @@ bandwidths(C::ConcreteConversion{<:Ultraspherical,<:Chebyshev}) =
bandwidths(C::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) =
0,order(domainspace(C))+1==order(rangespace(C)) ? 2 : ℵ₀

Base.stride(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical{Int}}) = 2
Base.stride(C::ConcreteConversion{<:Chebyshev,<:Ultraspherical{<:Union{Integer,StaticInt}}}) = 2
Base.stride(C::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) = 2

isdiag(::ConcreteConversion{<:Chebyshev,<:Ultraspherical}) = false
Expand All @@ -239,15 +243,20 @@ isdiag(::ConcreteConversion{<:Ultraspherical,<:Ultraspherical}) = false
## coefficients

# return the space that has banded Conversion to the other
function conversion_rule(a::Chebyshev,b::Ultraspherical{Int})
function conversion_rule(a::Chebyshev,b::Ultraspherical{<:Union{Integer,StaticInt}})
if domainscompatible(a,b)
a
else
NoSpace()
end
end

conversion_rule(a::Ultraspherical{<:StaticInt}, b::Ultraspherical{<:StaticInt}) =
_conversion_rule(a, b)
function conversion_rule(a::Ultraspherical{LT},b::Ultraspherical{LT}) where {LT}
_conversion_rule(a, b)
end
function _conversion_rule(a::Ultraspherical, b::Ultraspherical)
if domainscompatible(a,b) && isapproxinteger(order(a)-order(b))
order(a) < order(b) ? a : b
else
Expand All @@ -256,7 +265,7 @@ function conversion_rule(a::Ultraspherical{LT},b::Ultraspherical{LT}) where {LT}
end


function coefficients(g::AbstractVector,sp::Ultraspherical{Int},C::Chebyshev)
function coefficients(g::AbstractVector,sp::Ultraspherical{<:Union{Integer,StaticInt}},C::Chebyshev)
domainscompatible(C,sp) || throw(ArgumentError("domains don't match: $(domain(sp)) and $(domain(C))"))
if order(sp) == 1
ultraiconversion(g)
Expand Down Expand Up @@ -286,8 +295,7 @@ end

## Non-banded conversions

function getindex(M::ConcreteConversion{C,Ultraspherical{LT,DD,RR},T},
k::Integer,j::Integer) where {DD,RR,LT,C<:Chebyshev,T}
function getindex(M::ConcreteConversion{<:Chebyshev,<:Ultraspherical,T}, k::Integer,j::Integer) where {T}
λ = order(rangespace(M))
if λ == 1
if k==j==1
Expand Down Expand Up @@ -317,8 +325,7 @@ function getindex(M::ConcreteConversion{C,Ultraspherical{LT,DD,RR},T},
end


function getindex(M::ConcreteConversion{Ultraspherical{LT,DD,RR},C,T},
k::Integer,j::Integer) where {DD,RR,LT,C<:Chebyshev,T}
function getindex(M::ConcreteConversion{<:Ultraspherical,<:Chebyshev,T}, k::Integer,j::Integer) where {T}
λ = order(domainspace(M))
if λ == 1
# order must be 1
Expand Down
Loading