Support for cumsum with ChebyshevT (#53)

* Support for cumsum with ChebyshevT

* BigFloat with Ultraspherical

* Support ContinuumArrays v0.9

* increase coverage

Author: dlfivefifty

Project.toml

@@ -29,7 +29,7 @@ ArrayLayouts = "0.7"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.15, 0.16"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.8.5"
+ContinuumArrays = "0.8.5, 0.9"
 DomainSets = "0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"

examples/ultrasphericalspectralmethod.jl

@@ -1,4 +1,5 @@
 using ClassicalOrthogonalPolynomials, Plots, Test
+import ArrayLayouts: diagonal
 
 ###
 # We can solve ODEs like the Airy equation
@@ -14,8 +15,27 @@ C = Ultraspherical(2)
 x = axes(T,1)
 D = Derivative(x)
 
-c = [T[[begin,end],:]; C \ (D^2 * T - x .* T)] \ [airyai(-1); airyai(1); zeros(∞)]
+c = [T[[begin,end],:]; C \ ((D^2 - diagonal(x))*T)] \ [airyai(-1); airyai(1); zeros(∞)]
 u = T*c
 
 @test u[0.0] ≈ airyai(0.0)
-plot(u)
+plot(u)
+
+
+##
+# Lee & Greengard
+# ε*u'' - x*u' + u = 0, u(-1) = 1, u(1) = 2
+##
+
+T = ChebyshevT()
+C = Ultraspherical(2)
+x = axes(T,1)
+D = Derivative(x)
+
+ε = 1/100
+A = [T[[begin,end],:]; C \ ((ε*D^2 - x .* D + I) * T)]
+c = A \ [1; 2; zeros(∞)]
+u = T*c
+plot(u)
+
+C \ (ε*D^2 - x .* D + I) 

src/ClassicalOrthogonalPolynomials.jl

@@ -7,7 +7,7 @@ using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, Bloc
 
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
-                first, last, Slice, size, length, axes, IdentityUnitRange, sum, _sum,
+                first, last, Slice, size, length, axes, IdentityUnitRange, sum, _sum, cumsum,
                 to_indices, _maybetail, tail, getproperty, inv, show, isapprox, summary
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, adjointlayout,
@@ -19,7 +19,7 @@ import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_mat
 import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, AbstractLazyBandedLayout
 import LinearAlgebra: pinv, factorize, qr, adjoint, transpose, dot
 import BandedMatrices: AbstractBandedLayout, AbstractBandedMatrix, _BandedMatrix, bandeddata
-import FillArrays: AbstractFill, getindex_value
+import FillArrays: AbstractFill, getindex_value, SquareEye
 
 import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                     QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,

src/classical/chebyshev.jl

@@ -158,6 +158,8 @@ end
 # Conversion
 #####
 
+\(::Chebyshev{kind,T}, ::Chebyshev{kind,V}) where {kind,T,V} = SquareEye{promote_type(T,V)}(ℵ₀)
+
 function \(U::ChebyshevU, C::ChebyshevT)
     T = promote_type(eltype(U), eltype(C))
     _BandedMatrix(Vcat(-Ones{T}(1,∞)/2,
@@ -218,6 +220,8 @@ function \(A::ChebyshevT, B::Legendre)
         end, 1:∞, (1:∞)'))
 end
 
+\(A::AbstractJacobi, B::Chebyshev) = ApplyArray(inv,B \ A)
+
 
 function \(A::Jacobi, B::WeightedBasis{<:Any,<:JacobiWeight,<:Chebyshev})
     w, T = B.args
@@ -253,6 +257,13 @@ function _sum(A::WeightedBasis{T,<:ChebyshevWeight,<:Chebyshev}, dims) where T
     Hcat(convert(T, π), Zeros{T}(1,∞))
 end
 
+function cumsum(T::ChebyshevT{V}; dims::Integer) where V
+    @assert dims == 1
+    Σ = _BandedMatrix(Vcat(-one(V) ./ (-2:2:∞)', Zeros{V}(1,∞), Hcat(one(V), one(V) ./ (4:2:∞)')), ℵ₀, 0, 2)
+    ApplyQuasiArray(*, T, Vcat((-1).^(0:∞)'* Σ, Σ))
+end
+
+cumsum(f::Expansion{<:Any,<:ChebyshevT}) = cumsum(f.args[1]; dims=1) * f.args[2]
 
 ####
 # algebra

src/classical/jacobi.jl

@@ -516,4 +516,6 @@ function _sum(P::Legendre{T}, dims) where T
 end
 
 _sum(p::SubQuasiArray{T,1,Legendre{T},<:Tuple{Inclusion,Int}}, ::Colon) where T = parentindices(p)[2] == 1 ? convert(T, 2) : zero(T)
+_sum(P::AbstractJacobi{T}, dims) where T = 2 * (Legendre{T}() \ P)[1:1,:]
+
 

src/classical/ultraspherical.jl

@@ -73,8 +73,8 @@ function jacobimatrix(P::Ultraspherical{T}) where T
 end
 
 # These return vectors A[k], B[k], C[k] are from DLMF. Cause of MikaelSlevinsky we need an extra entry in C ... for now.
-function recurrencecoefficients(C::Ultraspherical)
-    λ = C.λ
+function recurrencecoefficients(C::Ultraspherical{T}) where T
+    λ = convert(T,C.λ)
     n = 0:∞
     (2(n .+ λ) ./ (n .+ 1), Zeros{typeof(λ)}(∞), (n .+ (2λ-1)) ./ (n .+ 1))
 end
@@ -91,28 +91,29 @@ end
 @simplify function *(D::Derivative{<:Any,<:ChebyshevInterval}, S::Legendre)
     T = promote_type(eltype(D),eltype(S))
     A = _BandedMatrix(Ones{T}(1,∞), ℵ₀, -1,1)
-    ApplyQuasiMatrix(*, Ultraspherical{T}(3/2), A)
+    ApplyQuasiMatrix(*, Ultraspherical{T}(convert(T,3)/2), A)
 end
 
 
 # Ultraspherical(λ+1)\(D*Ultraspherical(λ))
 @simplify function *(D::Derivative{<:Any,<:ChebyshevInterval}, S::Ultraspherical)
-    A = _BandedMatrix(Fill(2S.λ,1,∞), ℵ₀, -1,1)
-    ApplyQuasiMatrix(*, Ultraspherical{eltype(S)}(S.λ+1), A)
+    T = promote_type(eltype(D),eltype(S))
+    A = _BandedMatrix(Fill(2convert(T,S.λ),1,∞), ℵ₀, -1,1)
+    ApplyQuasiMatrix(*, Ultraspherical{T}(S.λ+1), A)
 end
 
 # Ultraspherical(λ-1)\ (D*wUltraspherical(λ))
 @simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::Weighted{<:Any,<:Ultraspherical})
     S = WS.P
     λ = S.λ
-    T = eltype(WS)
+    T = promote_type(eltype(D),eltype(WS))
     if λ == 1
-        A = _BandedMatrix((-(1:∞))', ℵ₀, 1,-1)
-        ApplyQuasiMatrix(*, ChebyshevTWeight{T}() .* ChebyshevT{T}(), A)
+        A = _BandedMatrix((-(one(T):∞))', ℵ₀, 1,-1)
+        ApplyQuasiMatrix(*, Weighted(ChebyshevT{T}()), A)
     else
         n = (0:∞)
         A = _BandedMatrix((-one(T)/(2*(λ-1)) * ((n.+1) .* (n .+ (2λ-1))))', ℵ₀, 1,-1)
-        ApplyQuasiMatrix(*, WeightedUltraspherical{T}(λ-1), A)
+        ApplyQuasiMatrix(*, Weighted(Ultraspherical{T}(λ-1)), A)
     end
 end
 
@@ -160,30 +161,32 @@ end
 \(T::Chebyshev, C::Ultraspherical) = inv(C \ T)
 
 function \(C2::Ultraspherical{<:Any,<:Integer}, C1::Ultraspherical{<:Any,<:Integer})
-    λ = C1.λ
     T = promote_type(eltype(C2), eltype(C1))
-    if C2.λ == λ+1
-        _BandedMatrix( Vcat(-(λ ./ ((0:∞) .+ λ))', Zeros(1,∞), (λ ./ ((0:∞) .+ λ))'), ℵ₀, 0, 2)
-    elseif C2.λ == λ
+    λ_Int = C1.λ
+    λ = convert(T,λ_Int)
+    if C2.λ == λ_Int+1
+        _BandedMatrix( Vcat(-(λ ./ ((0:∞) .+ λ))', Zeros{T}(1,∞), (λ ./ ((0:∞) .+ λ))'), ℵ₀, 0, 2)
+    elseif C2.λ == λ_Int
         Eye{T}(∞)
-    elseif C2.λ > λ
-        (C2 \ Ultraspherical(λ+1)) * (Ultraspherical(λ+1)\C1)
+    elseif C2.λ > λ_Int
+        (C2 \ Ultraspherical(λ_Int+1)) * (Ultraspherical(λ_Int+1)\C1)
     else
         error("Not implemented")
     end
 end
 
 function \(C2::Ultraspherical, C1::Ultraspherical)
-    λ = C1.λ
     T = promote_type(eltype(C2), eltype(C1))
+    λ_Int = C1.λ
+    λ = convert(T,λ_Int)
     if C2.λ == λ+1
-        _BandedMatrix( Vcat(-(λ ./ ((0:∞) .+ λ))', Zeros(1,∞), (λ ./ ((0:∞) .+ λ))'), ℵ₀, 0, 2)
-    elseif C2.λ == λ
+        _BandedMatrix( Vcat(-(λ ./ ((0:∞) .+ λ))', Zeros{T}(1,∞), (λ ./ ((0:∞) .+ λ))'), ℵ₀, 0, 2)
+    elseif C2.λ == λ_Int
         Eye{T}(∞)
-    elseif isinteger(C2.λ-λ) && C2.λ > λ
-        Cm = Ultraspherical{T}(λ+1)
+    elseif isinteger(C2.λ-λ_Int) && C2.λ > λ_Int
+        Cm = Ultraspherical{T}(λ_Int+1)
         (C2 \ Cm) * (Cm \ C1)
-    elseif isinteger(C2.λ-λ)
+    elseif isinteger(C2.λ-λ_Int)
         inv(C1 \ C2)
     else
         error("Not implemented")
@@ -193,13 +196,14 @@ end
 function \(w_A::WeightedUltraspherical, w_B::WeightedUltraspherical)
     wA,A = w_A.args
     wB,B = w_B.args
+    T = promote_type(eltype(w_A),eltype(w_B))
 
     if wA == wB
         A \ B
     elseif B.λ == A.λ+1 && wB.λ == wA.λ+1 # Lower
-        λ = A.λ
+        λ = convert(T,A.λ)
         _BandedMatrix(Vcat(((2λ:∞) .* ((2λ+1):∞) ./ (4λ .* (λ+1:∞)))',
-                            Zeros(1,∞),
+                            Zeros{T}(1,∞),
                             (-(1:∞) .* (2:∞) ./ (4λ .* (λ+1:∞)))'), ℵ₀, 2,0)
     else
         error("not implemented for $A and $wB")

src/lanczos.jl

@@ -232,6 +232,7 @@ struct LanczosLayout <: AbstractBasisLayout end
 MemoryLayout(::Type{<:LanczosPolynomial}) = LanczosLayout()
 arguments(::ApplyLayout{typeof(*)}, Q::LanczosPolynomial) = Q.P, LanczosConversion(Q.data)
 copy(L::Ldiv{LanczosLayout,Lay}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
+copy(L::Ldiv{LanczosLayout,Lay}) where Lay<:BroadcastLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
 copy(L::Ldiv{LanczosLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
 copy(L::Ldiv{LanczosLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
 LazyArrays._mul_arguments(Q::LanczosPolynomial) = arguments(ApplyLayout{typeof(*)}(), Q)

src/normalized.jl

@@ -182,6 +182,7 @@ function copy(L::Ldiv{WeightedOPLayout,WeightedOPLayout})
     convert(WeightedOrthogonalPolynomial, L.A) \ convert(WeightedOrthogonalPolynomial, L.B)
 end
 
+copy(L::Ldiv{WeightedOPLayout,<:BroadcastLayout}) = convert(WeightedOrthogonalPolynomial, L.A) \ L.B
 copy(L::Ldiv{WeightedOPLayout,<:AbstractLazyLayout}) = convert(WeightedOrthogonalPolynomial, L.A) \ L.B
 copy(L::Ldiv{WeightedOPLayout,<:AbstractBasisLayout}) = convert(WeightedOrthogonalPolynomial, L.A) \ L.B
 copy(L::Ldiv{<:AbstractLazyLayout,WeightedOPLayout}) = L.A \ convert(WeightedOrthogonalPolynomial, L.B)

test/test_chebyshev.jl

@@ -245,6 +245,18 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
         x = Inclusion(0..1)
         @test sum(wT[2x .- 1, :]; dims=1)[1,1:10] == [π/2; zeros(9)]
         @test sum(wT[2x .- 1, :] * [[1,2,3]; zeros(∞)]) == π/2
+
+        T = Chebyshev()
+        x = axes(T,1)
+        @test sum(T * (T \ exp.(x))) ≈ ℯ - inv(ℯ)
+    end
+
+    @testset "cumsum" begin
+        T = ChebyshevT()
+        x = axes(T,1)
+        @test (T \ cumsum(T; dims=1)) * (T \ exp.(x)) ≈ T \ (exp.(x) .- exp(-1))
+        f = T * (T \ exp.(x))
+        @test T \ cumsum(f) ≈ T \ (exp.(x) .- exp(-1))
     end
 
     @testset "algebra" begin

test/test_ultraspherical.jl

@@ -118,4 +118,17 @@ import LazyArrays: rowsupport, colsupport
         @test Ultraspherical(1/2) \ (JacobiWeight(0,0) .* Jacobi(0,0)) isa Diagonal
         @test (JacobiWeight(0,0) .* Jacobi(0,0)) \ Ultraspherical(1/2) isa Diagonal
     end
+
+    @testset "BigFloat" begin
+        U = Ultraspherical{BigFloat}(1)
+        T = ChebyshevT{BigFloat}()
+        x = axes(U,1)
+        D = Derivative(x)
+        
+        @test Weighted(T) \ (D * Weighted(U)) isa BandedMatrix{BigFloat}
+
+        C³ = Ultraspherical{BigFloat}(3)
+        c = [1; 2; 3; zeros(BigFloat,∞)]
+        @test C³[big(1)/10,:]'*(C³ \ U) * c ≈ U[big(1)/10,:]'c
+    end
 end