Add roots and extend Jacobi recurrences (#114)

* Add roots and extend Jacobi recurrences

* searchsortedfirst

* add tests

* checkpoints in Fourier should use axes type not eltype

* use axestype for grid

* Start Laurent

* Support Laurent

* Update Project.toml

* new grid gramework

* Update Project.toml

* Create annulipdes.jl

* Update annulipdes.jl

* arr -> szs

* tests pass

* Update ClassicalOrthogonalPolynomials.jl

* Update Project.toml

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.6.9"
+version = "0.7"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,8 +29,8 @@ ArrayLayouts = "0.8"
 BandedMatrices = "0.17"
 BlockArrays = "0.16.9"
 BlockBandedMatrices = "0.11.6"
-ContinuumArrays = "0.11"
-DomainSets = "0.5.6"
+ContinuumArrays = "0.12.1"
+DomainSets = "0.5.6, 0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3, 0.5"
 FastTransforms = "0.14.4"
@@ -39,8 +39,8 @@ HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.12.3"
 InfiniteLinearAlgebra = "0.6.7"
 IntervalSets = "0.5, 0.6, 0.7"
-LazyArrays = "0.22.11"
-LazyBandedMatrices = "0.7.14, 0.8"
+LazyArrays = "0.22.16"
+LazyBandedMatrices = "0.8.5"
 QuasiArrays = "0.9"
 SpecialFunctions = "1.0, 2"
 julia = "1.7"

examples/annulipdes.jl

@@ -0,0 +1,44 @@
+using ClassicalOrthogonalPolynomials, Plots
+
+ρ = 0.5; T = chebyshevt(ρ..1); U = chebyshevu(T); C = ultraspherical(2, ρ..1); r = axes(T,1); D = Derivative(r);
+
+L = C \ (r.^2 .* (D^2 * T)) + C \ (r .* (D * T)) # r^2 * ∂^2 + r*∂
+M = C\T # Identity
+
+m = 2
+Δₘ = L - m^2*M # r^2 * Laplacian for exp(im*m*θ)*u(r), i.e. (r^2 * ∂^2 + r*∂ - m^2*I)*u
+
+
+
+# Poisson solve
+c = C \ exp.(r)
+R = C \ (r .* C)
+
+
+d = [T[[begin,end],:];
+            Δₘ] \ [1; 2; R^2 * c]
+
+plot(T*d)
+
+
+# Helmholtz
+
+
+Q = R^2 * M # r^2, needed for Helmholtz (Δ + k^2)*u = f
+
+k = 5 # f
+
+d = [T[[begin,end],:];
+            Δₘ+k^2*Q] \ [1; 2; R^2 * c]
+
+
+# transform
+
+f = (r,θ) -> exp(r*cos(θ))
+
+n = 10 # create a 10 x 10 transform
+
+F = Fourier()
+𝐫,𝛉 = grid(T, n),grid(F, n)
+
+transform(T, transform(F, f.(𝐫, 𝛉'); dims=2); dims=1)

src/ClassicalOrthogonalPolynomials.jl

@@ -10,7 +10,8 @@ using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, Bloc
 import Base: @_inline_meta, axes, getindex, unsafe_getindex, convert, prod, *, /, \, +, -,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, setindex,
                 first, last, Slice, size, length, axes, IdentityUnitRange, sum, _sum, cumsum,
-                to_indices, _maybetail, tail, getproperty, inv, show, isapprox, summary
+                to_indices, _maybetail, tail, getproperty, inv, show, isapprox, summary,
+                findall, searchsortedfirst
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted, Broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, adjointlayout,
                 sub_materialize, arguments, sub_paddeddata, paddeddata, PaddedLayout, resizedata!, LazyVector, ApplyLayout, call,
@@ -34,10 +35,10 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import InfiniteLinearAlgebra: chop!, chop, choplength, compatible_resize!
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
-    inbounds_getindex, grid, plotgrid, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
+    inbounds_getindex, grid, plotgrid, _plotgrid, _grid, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
     AffineQuasiVector, AffineMap, WeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
     checkpoints, weight, unweighted, MappedBasisLayouts, __sum, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, _broadcastbasis,
-    plan_transform, plan_grid_transform
+    plan_transform, plan_grid_transform, MAX_PLOT_POINTS
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints, Plan
 
@@ -47,13 +48,14 @@ import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray, block
 import BandedMatrices: bandwidths
 
 export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomial,
-            Hermite, Jacobi, Legendre, Chebyshev, ChebyshevT, ChebyshevU, ChebyshevInterval, Ultraspherical, Fourier, Laguerre,
+            Hermite, Jacobi, Legendre, Chebyshev, ChebyshevT, ChebyshevU, ChebyshevInterval, Ultraspherical, Fourier, Laurent, Laguerre,
             HermiteWeight, JacobiWeight, ChebyshevWeight, ChebyshevGrid, ChebyshevTWeight, ChebyshevUWeight, UltrasphericalWeight, LegendreWeight, LaguerreWeight,
             WeightedUltraspherical, WeightedChebyshev, WeightedChebyshevT, WeightedChebyshevU, WeightedJacobi,
             ∞, Derivative, .., Inclusion,
             chebyshevt, chebyshevu, legendre, jacobi, ultraspherical,
             legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh, normalizedjacobip,
-            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted, PiecewiseInterlace, plan_transform
+            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted, PiecewiseInterlace, plan_transform,
+            expand, transform
 
 
 import Base: oneto
@@ -111,7 +113,6 @@ copy(L::Ldiv{MappedOPLayout,Lay}) where Lay<:MappedBasisLayouts = copy(Ldiv{Mapp
 
 # OPs are immutable
 copy(a::OrthogonalPolynomial) = a
-copy(a::SubQuasiArray{<:Any,N,<:OrthogonalPolynomial}) where N = a
 
 """
     jacobimatrix(P)
@@ -245,18 +246,10 @@ _tritrunc(X, n) = _tritrunc(MemoryLayout(X), X, n)
 jacobimatrix(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Inclusion,OneTo}}) =
     _tritrunc(jacobimatrix(parent(V)), maximum(parentindices(V)[2]))
 
-grid(P::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,OneTo}}) =
-    eigvals(symtridiagonalize(jacobimatrix(P)))
-
-function grid(Pn::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,Any}})
-    kr,jr = parentindices(Pn)
-    grid(parent(Pn)[:,oneto(maximum(jr))])
-end
-
-function plotgrid(Pn::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,Any}}) where T
-    kr,jr = parentindices(Pn)
-    grid(parent(Pn)[:,oneto(40maximum(jr))])
-end
+_grid(::AbstractOPLayout, P, n::Integer) = eigvals(symtridiagonalize(jacobimatrix(P[:,OneTo(n)])))
+_grid(::MappedOPLayout, P, n::Integer) = _grid(MappedBasisLayout(), P, n)
+_plotgrid(::AbstractOPLayout, P, n::Integer) = grid(P, min(40n, MAX_PLOT_POINTS))
+_plotgrid(::MappedOPLayout, P, n::Integer) = _plotgrid(MappedBasisLayout(), P, n)
 
 function golubwelsch(X)
     D, V = eigen(symtridiagonalize(X))  # Eigenvalue decomposition
@@ -324,16 +317,16 @@ end
 *(A::AbstractMatrix, P::MulPlan) = MulPlan(A*P.matrix, P.dims)
 
 
-function plan_grid_transform(Q::Normalized, arr, dims=1:ndims(arr))
-    L = Q[:,OneTo(size(arr,1))]
+function plan_grid_transform(Q::Normalized, szs::NTuple{N,Int}, dims=1:N) where N
+    L = Q[:,OneTo(szs[1])]
     x,w = golubwelsch(L)
     x, MulPlan(L[x,:]'*Diagonal(w), dims)
 end
 
-function plan_grid_transform(P::OrthogonalPolynomial, arr, dims...)
+function plan_grid_transform(P::OrthogonalPolynomial, szs::NTuple{N,Int}, dims=1:N) where N
     Q = Normalized(P)
-    x, A = plan_grid_transform(Q, arr, dims...)
-    n = size(arr,1)
+    x, A = plan_grid_transform(Q, szs, dims...)
+    n = szs[1]
     D = (P \ Q)[1:n, 1:n]
     x, D * A
 end
@@ -381,6 +374,6 @@ include("classical/ultraspherical.jl")
 include("classical/laguerre.jl")
 include("classical/fourier.jl")
 include("stieltjes.jl")
-
+include("roots.jl")
 
 end # module

src/classical/chebyshev.jl

@@ -103,17 +103,17 @@ Jacobi(C::ChebyshevU{T}) where T = Jacobi(one(T)/2,one(T)/2)
 #######
 
 
-function plan_grid_transform(T::ChebyshevT, arr, dims...)
-    n = size(arr,1)
-    x = grid(T[:,oneto(n)])
+function plan_grid_transform(T::ChebyshevT, arr::AbstractArray, dims...)
+    x = grid(T, size(arr,1))
     x, plan_chebyshevtransform(arr, dims...)
 end
-function plan_grid_transform(U::ChebyshevU{<:FastTransforms.fftwNumber}, arr, dims...)
-    n = size(arr,1)
-    x = grid(U[:,oneto(n)])
+function plan_grid_transform(U::ChebyshevU{<:FastTransforms.fftwNumber}, arr::AbstractArray, dims...)
+    x = grid(U, size(arr,1))
     x, plan_chebyshevutransform(arr, dims...)
 end
 
+plan_grid_transform(T::Chebyshev, szs::NTuple{N,Int}, dims...) where N = plan_grid_transform(T, Array{eltype(T)}(undef, szs...), dims...)
+
 
 ########
 # Jacobi Matrix

src/classical/fourier.jl

@@ -1,50 +1,69 @@
+abstract type AbstractFourier{T} <: Basis{T} end
+
+struct Fourier{T} <: AbstractFourier{T} end
+struct Laurent{T} <: AbstractFourier{T} end
 
-struct Fourier{T} <: Basis{T} end
 Fourier() = Fourier{Float64}()
+Laurent() = Laurent{ComplexF64}()
 
 ==(::Fourier, ::Fourier) = true
+==(::Laurent, ::Laurent) = true
 
-axes(F::Fourier) = (Inclusion(ℝ), _BlockedUnitRange(1:2:∞))
+axes(F::AbstractFourier) = (Inclusion(ℝ), _BlockedUnitRange(1:2:∞))
 
 function getindex(F::Fourier{T}, x::Real, j::Int)::T where T
     isodd(j) && return cos((j÷2)*x)
     sin((j÷2)*x)
 end
 
+function getindex(F::Laurent{T}, x::Real, j::Int)::T where T
+    s = 1-2iseven(j)
+    exp(im*s*(j÷2)*x)
+end
+
 ### transform
-checkpoints(::Fourier{T}) where T = T[1.223972,3.14,5.83273484]
+checkpoints(F::AbstractFourier) = eltype(axes(F,1))[1.223972,3.14,5.83273484]
 
 fouriergrid(T, n) = convert(T,π)*collect(0:2:2n-2)/n
+grid(Pn::AbstractFourier, n::Integer) = fouriergrid(eltype(axes(Pn,1)), n)
 
-function grid(Pn::SubQuasiArray{T,2,<:Fourier,<:Tuple{<:Inclusion,<:AbstractUnitRange}}) where T
-    kr,jr = parentindices(Pn)
-    n = maximum(jr)
-    fouriergrid(T, n)
-end
+abstract type AbstractShuffledPlan{T} <: Plan{T} end
 
 """
 Gives a shuffled version of the real FFT, with order
 1,sin(θ),cos(θ),sin(2θ)…
 """
-struct ShuffledRFFT{T,Pl<:Plan} <: Factorization{T}
+struct ShuffledR2HC{T,Pl<:Plan} <: AbstractShuffledPlan{T}
+    plan::Pl
+end
+
+"""
+Gives a shuffled version of the FFT, with order
+1,sin(θ),cos(θ),sin(2θ)…
+"""
+struct ShuffledFFT{T,Pl<:Plan} <: AbstractShuffledPlan{T}
     plan::Pl
 end
 
-size(F::ShuffledRFFT, k) = size(F.plan,k)
-size(F::ShuffledRFFT) = size(F.plan)
+size(F::AbstractShuffledPlan, k) = size(F.plan,k)
+size(F::AbstractShuffledPlan) = size(F.plan)
 
-ShuffledRFFT{T}(p::Pl) where {T,Pl<:Plan} = ShuffledRFFT{T,Pl}(p)
-ShuffledRFFT{T}(n, d...) where T = ShuffledRFFT{T}(FFTW.plan_r2r(Array{T}(undef, n), FFTW.R2HC, d...))
+ShuffledR2HC{T}(p::Pl) where {T,Pl<:Plan} = ShuffledR2HC{T,Pl}(p)
+ShuffledR2HC{T}(n, d...) where T = ShuffledR2HC{T}(FFTW.plan_r2r(Array{T}(undef, n), FFTW.R2HC, d...))
 
-function _shuffledrfft_postscale!(_, ret::AbstractVector{T}) where T
+ShuffledFFT{T}(p::Pl) where {T,Pl<:Plan} = ShuffledFFT{T,Pl}(p)
+ShuffledFFT{T}(n, d...) where T = ShuffledFFT{T}(FFTW.plan_fft(Array{T}(undef, n), d...))
+
+
+function _shuffledR2HC_postscale!(_, ret::AbstractVector{T}) where T
     n = length(ret)
     lmul!(convert(T,2)/n, ret)
     ret[1] /= 2
     iseven(n) && (ret[n÷2+1] /= 2)
     negateeven!(reverseeven!(interlace!(ret,1)))
 end
 
-function _shuffledrfft_postscale!(d::Number, ret::AbstractMatrix{T}) where T
+function _shuffledR2HC_postscale!(d::Number, ret::AbstractMatrix{T}) where T
     if isone(d)
         n = size(ret,1)
         lmul!(convert(T,2)/n, ret)
@@ -65,20 +84,53 @@ function _shuffledrfft_postscale!(d::Number, ret::AbstractMatrix{T}) where T
     ret
 end
 
+function _shuffledFFT_postscale!(_, ret::AbstractVector{T}) where T
+    n = length(ret)
+    cfs = lmul!(inv(convert(T,n)), ret)
+    reverseeven!(interlace!(cfs,1))
+end
 
-function mul!(ret::AbstractArray{T}, F::ShuffledRFFT{T}, b::AbstractArray) where T
+function _shuffledFFT_postscale!(d::Number, ret::AbstractMatrix{T}) where T
+    if isone(d)
+        n = size(ret,1)
+        lmul!(inv(convert(T,n)), ret)
+        for j in axes(ret,2)
+            reverseeven!(interlace!(view(ret,:,j),1))
+        end
+    else
+        n = size(ret,2)
+        lmul!(inv(convert(T,n)), ret)
+        for k in axes(ret,1)
+            reverseeven!(interlace!(view(ret,k,:),1))
+        end
+    end
+    ret
+end
+
+function mul!(ret::AbstractArray{T}, F::ShuffledR2HC{T}, b::AbstractArray) where T
+    mul!(ret, F.plan, convert(Array{T}, b))
+    _shuffledR2HC_postscale!(F.plan.region, ret)
+end
+
+function mul!(ret::AbstractArray{T}, F::ShuffledFFT{T}, b::AbstractArray) where T
     mul!(ret, F.plan, convert(Array{T}, b))
-    _shuffledrfft_postscale!(F.plan.region, ret)
+    _shuffledFFT_postscale!(F.plan.region, ret)
 end
 
-*(F::ShuffledRFFT{T}, b::AbstractVecOrMat) where T = mul!(similar(b, T), F, b)
+*(F::AbstractShuffledPlan{T}, b::AbstractVecOrMat) where T = mul!(similar(b, T), F, b)
 
 factorize(L::SubQuasiArray{T,2,<:Fourier,<:Tuple{<:Inclusion,<:OneTo}}) where T =
-    TransformFactorization(grid(L), ShuffledRFFT{T}(size(L,2)))
+    TransformFactorization(grid(L), ShuffledR2HC{T}(size(L,2)))
 factorize(L::SubQuasiArray{T,2,<:Fourier,<:Tuple{<:Inclusion,<:OneTo}}, d) where T =
-    TransformFactorization(grid(L), ShuffledRFFT{T}((size(L,2),d),1))
+    TransformFactorization(grid(L), ShuffledR2HC{T}((size(L,2),d),1))
+
+factorize(L::SubQuasiArray{T,2,<:Laurent,<:Tuple{<:Inclusion,<:OneTo}}) where T =
+    TransformFactorization(grid(L), ShuffledFFT{T}(size(L,2)))
+factorize(L::SubQuasiArray{T,2,<:Laurent,<:Tuple{<:Inclusion,<:OneTo}}, d) where T =
+    TransformFactorization(grid(L), ShuffledFFT{T}((size(L,2),d),1))
+
 
-factorize(L::SubQuasiArray{T,2,<:Fourier,<:Tuple{<:Inclusion,<:BlockSlice}},d...) where T =
+factorize(L::SubQuasiArray{T,2,<:AbstractFourier,<:Tuple{<:Inclusion,<:BlockSlice}},d...) where T =
     ProjectionFactorization(factorize(parent(L)[:,OneTo(size(L,2))],d...),parentindices(L)[2])
 
 import BlockBandedMatrices: _BlockSkylineMatrix
@@ -88,11 +140,21 @@ import BlockBandedMatrices: _BlockSkylineMatrix
     PseudoBlockArray(Diagonal(Vcat(2convert(TV,π),Fill(convert(TV,π),∞))), (axes(A,1),axes(B,2)))
 end
 
+@simplify function *(A::QuasiAdjoint{<:Any,<:Laurent}, B::Laurent)
+    TV = promote_type(eltype(A),eltype(B))
+    Diagonal(Fill(2convert(TV,π),(axes(B,2),)))
+end
+
 @simplify function *(D::Derivative, F::Fourier)
     TV = promote_type(eltype(D),eltype(F))
     Fourier{TV}()*_BlockArray(Diagonal(Vcat([reshape([0.0],1,1)], (1.0:∞) .* Fill([0 -one(TV); one(TV) 0], ∞))), (axes(F,2),axes(F,2)))
 end
 
+@simplify function *(D::Derivative, F::Laurent)
+    TV = promote_type(eltype(D),eltype(F))
+    Laurent{TV}() * Diagonal(PseudoBlockVector((((1:∞) .÷ 2) .* (1 .- 2 .* iseven.(1:∞))) * convert(TV,im), (axes(F,2),)))
+end
+
 
 function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), c::BroadcastQuasiVector{<:Any,typeof(cos),<:Tuple{<:Inclusion{<:Any,RealNumbers}}}, F::Fourier)
     axes(c,1) == axes(F,1) || throw(DimensionMismatch())