up to ODETest

Author: dlfivefifty

src/ApproxFunOrthogonalPolynomials.jl

@@ -1,6 +1,6 @@
 module ApproxFunOrthogonalPolynomials
-using Base, LinearAlgebra, Reexport, AbstractFFTs, FFTW, InfiniteArrays, FillArrays, FastTransforms, IntervalSets, 
-            DomainSets, Statistics
+using Base, LinearAlgebra, Reexport, BandedMatrices, AbstractFFTs, FFTW, InfiniteArrays, FillArrays, FastTransforms, IntervalSets, 
+            DomainSets, Statistics, SpecialFunctions
 
 @reexport using ApproxFunBase
 
@@ -12,18 +12,20 @@ import ApproxFunBase: normalize!, flipsign, FiniteRange, Fun, MatrixFun, UnsetSp
                     UnivariateSpace, AmbiguousSpace, SumSpace, SubSpace, WeightSpace, NoSpace, Space,
                     HeavisideSpace, PointSpace,
                     IntervalOrSegment, RaggedMatrix, AlmostBandedMatrix,
-                    AnyDomain, ZeroSpace, TrivialInterlacer, BlockInterlacer, 
+                    AnyDomain, ZeroSpace, ArraySpace, TrivialInterlacer, BlockInterlacer, 
                     AbstractTransformPlan, TransformPlan, ITransformPlan,
                     ConcreteConversion, ConcreteMultiplication, ConcreteDerivative, ConcreteIntegral,
                     ConcreteVolterra, Volterra, VolterraWrapper,
-                    MultiplicationWrapper, ConversionWrapper, DerivativeWrapper, Evaluation,
+                    MultiplicationWrapper, ConversionWrapper, DerivativeWrapper, Evaluation, EvaluationWrapper,
                     Conversion, Multiplication, Derivative, Integral, bandwidths, 
                     ConcreteEvaluation, ConcreteDefiniteLineIntegral, ConcreteDefiniteIntegral, ConcreteIntegral,
                     DefiniteLineIntegral, DefiniteIntegral, ConcreteDefiniteIntegral, ConcreteDefiniteLineIntegral,
                     ReverseOrientation, ReverseOrientationWrapper, ReverseWrapper, Reverse, NegateEven, Dirichlet, ConcreteDirichlet,
                     TridiagonalOperator, SubOperator, Space, @containsconstants, spacescompatible,
                     hasfasttransform, canonicalspace, domain, setdomain, prectype, domainscompatible, 
-                    plan_transform, plan_itransform, plan_transform!, plan_itransform!, transform, itransform, hasfasttransform, Integral, 
+                    plan_transform, plan_itransform, plan_transform!, plan_itransform!, transform, itransform, hasfasttransform, 
+                    CanonicalTransformPlan,
+                    Integral, 
                     domainspace, rangespace, boundary, 
                     union_rule, conversion_rule, maxspace_rule, conversion_type, maxspace, hasconversion, points, 
                     rdirichlet, ldirichlet, lneumann, rneumann, ivp, bvp, 
@@ -34,14 +36,20 @@ import ApproxFunBase: normalize!, flipsign, FiniteRange, Fun, MatrixFun, UnsetSp
                     invfromcanonicalD, fromcanonical, tocanonical, fromcanonicalD, tocanonicalD, canonicaldomain, setcanonicaldomain, mappoint,
                     reverseorientation, checkpoints, evaluate, mul_coefficients, coefficients, isconvertible,
                     clenshaw, ClenshawPlan, sineshaw,
-                    toeplitz_getindex, toeplitz_axpy!, ToeplitzOperator, hankel_getindex, 
+                    toeplitz_getindex, toeplitz_axpy!, sym_toeplitz_axpy!, hankel_axpy!, ToeplitzOperator, SymToeplitzOperator, hankel_getindex, 
                     SpaceOperator, ZeroOperator, InterlaceOperator,
                     interlace!, reverseeven!, negateeven!, cfstype, pad!,
-                    extremal_args, hesseneigvals, chebyshev_clenshaw, recA, recB, recC, roots, chebmult_getindex, intpow, alternatingsum
+                    extremal_args, hesseneigvals, chebyshev_clenshaw, recA, recB, recC, roots, chebmult_getindex, intpow, alternatingsum,
+                    domaintype, diagindshift, rangetype, weight, isapproxinteger, default_Dirichlet, scal!
 
 import DomainSets: Domain, indomain, UnionDomain, ProductDomain, FullSpace, Point, elements, DifferenceDomain,
             Interval, ChebyshevInterval, boundary, ∂, rightendpoint, leftendpoint,
-            dimension, Domain1d, Domain2d         
+            dimension, Domain1d, Domain2d   
+            
+import BandedMatrices: bandrange, bandshift,
+                inbands_getindex, inbands_setindex!, bandwidth, AbstractBandedMatrix,
+                colstart, colstop, colrange, rowstart, rowstop, rowrange,
+                bandwidths, _BandedMatrix, BandedMatrix            
 
 import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !, !=, eltype, iterate,
                 >=, /, ^, \, ∪, transpose, size, reindex, tail, broadcast, broadcast!, copyto!, copy, to_index, (:),
@@ -66,10 +74,26 @@ import FastTransforms: ChebyshevTransformPlan, IChebyshevTransformPlan, plan_che
                         plan_chebyshevtransform!, plan_ichebyshevtransform, plan_ichebyshevtransform!,
                         pochhammer
 
+# we need to import all special functions to use Calculus.symbolic_derivatives_1arg
+# we can't do importall Base as we replace some Base definitions
+import SpecialFunctions: sinpi, cospi, airy, besselh,
+                    asinh, acosh,atanh, erfcx, dawson, erf, erfi,
+                    sin, cos, sinh, cosh, airyai, airybi, airyaiprime, airybiprime,
+                    hankelh1, hankelh2, besselj, besselj0, bessely, besseli, besselk,
+                    besselkx, hankelh1x, hankelh2x, exp2, exp10, log2, log10,
+                    tan, tanh, csc, asin, acsc, sec, acos, asec,
+                    cot, atan, acot, sinh, csch, asinh, acsch,
+                    sech, acosh, asech, tanh, coth, atanh, acoth,
+                    expm1, log1p, lfact, sinc, cosc, erfinv, erfcinv, beta, lbeta,
+                    eta, zeta, gamma,  lgamma, polygamma, invdigamma, digamma, trigamma,
+                    abs, sign, log, expm1, tan, abs2, sqrt, angle, max, min, cbrt, log,
+                    atan, acos, asin, erfc, inv
+
 include("ultraspherical.jl")
 include("Domains/Domains.jl")
 include("Spaces/Spaces.jl")
 include("roots.jl")
 include("specialfunctions.jl")
+include("fastops.jl")
 
 end

src/Spaces/Chebyshev/ChebyshevOperators.jl

@@ -1,5 +1,5 @@
 
-recA(::Type{T},::Chebyshev,k) where {T} = 2one(T)
+recA(::Type{T},::Chebyshev,k) where {T} = k == 0 ? one(T) : 2one(T)
 recB(::Type{T},::Chebyshev,_) where {T} = zero(T)
 recC(::Type{T},::Chebyshev,k) where {T} = one(T)   # one(T) ensures we get correct type
 
@@ -38,7 +38,7 @@ function evaluatechebyshev(n::Integer,x::T) where T<:Number
 end
 
 
-function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},typeof(leftendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
+function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
     T=eltype(op)
     if op.order == 0
         ifelse(isodd(j),  # right rule
@@ -50,7 +50,7 @@ function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},typeof(leftendpoint)},
     end
 end
 
-function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},typeof(rightendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
+function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
     T=eltype(op)
     if op.order == 0
         one(T)
@@ -60,7 +60,7 @@ function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},typeof(rightendpoint)}
     end
 end
 
-function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},typeof(leftendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
+function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
     T=eltype(op)
     x = op.x
     d = domain(op)
@@ -86,7 +86,7 @@ function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},typeof(leftendpoint)},
 end
 
 
-function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},typeof(rightendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
+function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
     T=eltype(op)
     x = op.x
     d = domain(op)

src/Spaces/PolynomialSpace.jl

@@ -201,7 +201,7 @@ function BandedMatrix(S::SubOperator{T,ConcreteMultiplication{C,PS,T},
         # for e_1^⊤ U_x\a == a[1]*I-(α-J)*a[2]/β == (a[1]-α*a[2]/β)*I + J*a[2]/β
         # implying
         α,β=recα(T,sp,1),recβ(T,sp,1)
-        ret=Operator{T}(ApproxFun.Recurrence(M.space))[kr,jr]::BandedMatrix{T}
+        ret=Operator{T}(Recurrence(M.space))[kr,jr]::BandedMatrix{T}
         lmul!(a[2]/β,ret)
         shft=kr[1]-jr[1]
         ret[band(shft)] .+= a[1]-α*a[2]/β
@@ -211,7 +211,7 @@ function BandedMatrix(S::SubOperator{T,ConcreteMultiplication{C,PS,T},
     jkr=max(1,min(jr[1],kr[1])-(n-1)÷2):max(jr[end],kr[end])+(n-1)÷2
 
     #Multiplication is transpose
-    J=Operator{T}(ApproxFun.Recurrence(M.space))[jkr,jkr]
+    J=Operator{T}(Recurrence(M.space))[jkr,jkr]
 
     B=n-1  # final bandwidth
 
@@ -318,14 +318,14 @@ end
 function getindex(op::ConcreteEvaluation{J,typeof(leftendpoint)},kr::AbstractRange) where J<:PolynomialSpace
     sp=op.space
     T=eltype(op)
-
+    @assert op.order == 0
     forwardrecurrence(T,sp,kr.-1,-one(T))
 end
 
 function getindex(op::ConcreteEvaluation{J,typeof(rightendpoint)},kr::AbstractRange) where J<:PolynomialSpace
     sp=op.space
     T=eltype(op)
-
+    @assert op.order == 0
     forwardrecurrence(T,sp,kr.-1,one(T))
 end
 
@@ -334,8 +334,8 @@ function getindex(op::ConcreteEvaluation{J,TT},kr::AbstractRange) where {J<:Poly
     sp=op.space
     T=eltype(op)
     x=op.x
-
-    forwardrecurrence(T,sp,kr.-1,tocanonical(sp,x))
+    @assert op.order == 0
+    forwardrecurrence(T,sp,kr .- 1,tocanonical(sp,x))
 end
 
 

src/fastops.jl

@@ -0,0 +1,108 @@
+###
+# This file contains BLAS/BandedMatrix overrides for operators
+# that depend on the structure of BandedMatrix
+####
+
+
+
+
+#####
+# Conversions
+#####
+
+function BandedMatrix(S::SubOperator{T,ConcreteConversion{Chebyshev{DD,RR},Ultraspherical{Int,DD,RR},T},
+                              Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,DD,RR}
+    # we can assume order is 1
+    ret = BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))
+    kr,jr = parentindices(S)
+    dg = diagindshift(S)
+
+    @assert -bandwidth(ret,1) ≤ dg ≤ bandwidth(ret,2)-2
+
+    ret[band(dg)] .= 0.5
+    ret[band(dg+1)] .= 0.0
+    ret[band(dg+2)] .= -0.5
+
+    # correct first entry
+    if 1 in kr && 1 in jr
+        ret[1,1] = 1.0
+    end
+
+    ret
+end
+
+function BandedMatrix(V::SubOperator{T,ConcreteConversion{Ultraspherical{LT,DD,RR},Ultraspherical{LT,DD,RR},T},
+                                                                  Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,LT,DD,RR}
+
+    n,m = size(V)
+    V_l, V_u = bandwidths(V)
+    ret = BandedMatrix{eltype(V)}(undef, (n,m), (V_l,V_u))
+    kr,jr = parentindices(V)
+    dg = diagindshift(V)
+
+
+    λ = order(rangespace(parent(V)))
+    c = λ-one(T)
+
+    # need to drop columns
+
+
+
+    1-n ≤ dg ≤ m-1 && (ret[band(dg)] .= c./(jr[max(0,dg)+1:min(n+dg,m)] .- 2 .+ λ))
+    1-n ≤ dg+1 ≤ m-1 && (ret[band(dg+1)] .= 0)
+    1-n ≤ dg+2 ≤ m-1 && (ret[band(dg+2)] .= c./(2 .- λ .- jr[max(0,dg+2)+1:min(n+dg+2,m)]))
+
+    ret
+end
+
+
+#####
+# Derivatives
+#####
+
+
+
+function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Chebyshev{DD,RR},K,T},
+                                                     Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,K,DD,RR}
+
+    n,m = size(S)
+    ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
+    kr,jr = parentindices(S)
+    dg = diagindshift(S)
+
+    D = parent(S)
+    k = D.order
+    d = domain(D)
+
+    C=convert(T,pochhammer(one(T),k-1)/2*(4/(complexlength(d)))^k)
+
+
+    # need to drop columns
+
+
+    if 1-n ≤ dg+k ≤ m-1
+        ret[band(dg+k)] .= C.*(jr[max(0,dg+k)+1:min(n+dg+k,m)] .- one(T))
+    end
+
+    ret
+end
+
+
+function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Ultraspherical{LT,DD,RR},K,T},
+                                                  Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,K,DD,RR,LT}
+    n,m = size(S)
+    ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
+    kr,jr = parentindices(S)
+    dg = diagindshift(S)
+
+    D = parent(S)
+    k = D.order
+    λ = order(domainspace(D))
+    d = domain(D)
+
+    C = convert(T,pochhammer(one(T)*λ,k)*(4/(complexlength(d)))^k)
+    ret[band(dg+k)] .= C
+
+    ret
+end
+

test/ChebyshevTest.jl

@@ -1,7 +1,12 @@
 using ApproxFunOrthogonalPolynomials, ApproxFunBase, LinearAlgebra, Test
-    import ApproxFunBase: testspace
+    import ApproxFunBase: testspace, recA, recB, recC
+    import ApproxFunOrthogonalPolynomials: forwardrecurrence
 
 @testset "Chebyshev" begin
+    @testset "Forward recurrence" begin
+        @test forwardrecurrence(Float64,Chebyshev(),0:9,1.0) == ones(10)
+    end
+
     @testset "ChebyshevInterval" begin
         @test Fun(x->2,10)(.1) ≈ 2
         @test Fun(x->2)(.1) ≈ 2
@@ -151,6 +156,7 @@ using ApproxFunOrthogonalPolynomials, ApproxFunBase, LinearAlgebra, Test
         x = Fun(identity,0..10)
         f = sin(x^2)
         g = cos(x)
+        @test exp(x)(0.1) ≈ exp(0.1)
         @test f(.1) ≈ sin(.1^2)
 
         x = Fun(identity,0..100)
@@ -160,7 +166,7 @@ using ApproxFunOrthogonalPolynomials, ApproxFunBase, LinearAlgebra, Test
 
     @testset "Reverse" begin
         f=Fun(exp)
-        @test ApproxFunOrthogonalPolynomials.default_Fun(f, Chebyshev(Segment(1 , -1)), ncoefficients(f))(0.1) ≈ exp(0.1)
+        @test ApproxFunBase.default_Fun(f, Chebyshev(Segment(1 , -1)), ncoefficients(f))(0.1) ≈ exp(0.1)
         @test Fun(f,Chebyshev(Segment(1,-1)))(0.1) ≈ f(0.1)
     end
 

test/ComplexTest.jl

@@ -1,4 +1,4 @@
-using ApproxFun, Test
+using ApproxFunOrthogonalPolynomials, Test
 
 @testset "Complex" begin
     @testset "Differentiation" begin

test/EigTest.jl

@@ -1,5 +1,5 @@
-using ApproxFun, BandedMatrices, SpecialFunctions, Test
-    import ApproxFun: bandwidth
+using ApproxFunOrthogonalPolynomials, ApproxFunBase, BandedMatrices, SpecialFunctions, Test
+    import ApproxFunBase: bandwidth
 
 #
 # Fix for BigFloat stackoverflow