start Jacobi

Author: dlfivefifty

REQUIRE

@@ -1,6 +1,7 @@
 julia 0.7
 IntervalSets 0.3.1
 DomainSets 0.0.1
+FillArrays 0.3
 LazyArrays 0.4
 BandedMatrices 0.7.2
 InfiniteArrays 0.0.2

src/ContinuumArrays.jl

@@ -1,16 +1,17 @@
 module ContinuumArrays
 using IntervalSets, LinearAlgebra, LazyArrays, BandedMatrices, InfiniteArrays, DomainSets
 import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
-                IndexStyle, IndexLinear, ==
+                IndexStyle, IndexLinear, ==, OneTo
 import Base.Broadcast: materialize
 import LazyArrays: Mul2
-import BandedMatrices: AbstractBandedLayout
+import BandedMatrices: AbstractBandedLayout, _BandedMatrix
 
 include("QuasiArrays/QuasiArrays.jl")
 using .QuasiArrays
-import .QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, slice, QSlice, SubQuasiArray
+import .QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, slice, QSlice, SubQuasiArray,
+                    QuasiDiagonal
 
-export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative
+export Spline, LinearSpline, HeavisideSpline, DiracDelta, Derivative, JacobiWeight, Jacobi, Legendre
 
 ####
 # Interval indexing support
@@ -35,6 +36,17 @@ function checkindex(::Type{Bool}, inds::AbstractInterval, I::AbstractArray)
     b
 end
 
+
+# we represent as a Mul with a banded matrix
+function materialize(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:QSlice,<:AbstractUnitRange}})
+    A = parent(V)
+    _,jr = parentindices(V)
+    first(jr) ≥ 1 || throw(BoundsError())
+    P = BandedMatrix((1-first(jr) => Ones{Int}(length(jr)),), (size(A,2), length(jr)))
+    A*P
+end
+
+
 include("operators.jl")
 include("bases/bases.jl")
 

src/QuasiArrays/matmul.jl

@@ -32,7 +32,7 @@ end
 QuasiMatMulMat{styleA, styleB, T, V} = QuasiArrayMulArray{styleA, styleB, 2, 2, T, V}
 QuasiMatMulQuasiMat{styleA, styleB, T, V} = QuasiArrayMulQuasiArray{styleA, styleB, 2, 2, T, V}
 
-
+*(A::AbstractQuasiArray, B::AbstractQuasiArray, C::AbstractQuasiArray) = materialize(Mul(A,B,C))
 *(A::AbstractQuasiArray, B::AbstractQuasiArray) = materialize(Mul(A,B))
 *(A::AbstractQuasiArray, B::AbstractArray) = materialize(Mul(A,B))
 *(A::AbstractArray, B::AbstractQuasiArray) = materialize(Mul(A,B))

src/QuasiArrays/quasidiagonal.jl

@@ -6,7 +6,6 @@ struct QuasiDiagonal{T,V<:AbstractQuasiVector{T}} <: AbstractQuasiMatrix{T}
     diag::V
 
     function QuasiDiagonal{T,V}(diag) where {T,V<:AbstractQuasiVector{T}}
-        @assert !has_offset_axes(diag)
         new{T,V}(diag)
     end
 end

src/bases/jacobi.jl

@@ -9,11 +9,35 @@ function getindex(w::JacobiWeight, x::Real)
     (1-x)^w.a * (1+x)^w.b
 end
 
+abstract type AbstractJacobi{T} <: AbstractQuasiMatrix{T} end
 
-struct Jacobi{T} <: AbstractQuasiMatrix{T}
+struct Legendre{T} <: AbstractJacobi{T} end
+Legendre() = Legendre{Float64}()
+
+struct Jacobi{T} <: AbstractJacobi{T}
     b::T
     a::T
 end
 
-axes(::Jacobi) = (ChebyshevInterval(), OneTo(∞))
+axes(::AbstractJacobi) = (ChebyshevInterval(), OneTo(∞))
 ==(P::Jacobi, Q::Jacobi) = P.a == Q.a && P.b == Q.b
+
+
+materialize(M::Mul2{<:Any,<:Any,<:QuasiAdjoint{<:Any,<:Legendre},<:Legendre}) =
+    Diagonal(2 ./ (2(0:∞) .+ 1))
+
+
+function materialize(M::Mul2{<:Any,<:Any,<:Derivative{<:Any,<:ChebyshevInterval},<:Jacobi})
+    _, S = M.factors
+    D = _BandedMatrix(((1:∞ .+ 1 .+ S.a .+ S.b)/2)', ∞, -1,1)
+    Mul(Jacobi(S.a+1,S.b+1), D)
+end
+
+function materialize(M::Mul{<:Tuple,<:Tuple{<:Derivative{<:Any,<:ChebyshevInterval},
+                                        QuasiDiagonal{Bool,JacobiWeight{Bool}},Jacobi{Bool}}})
+    _, W, S = M.factors
+    w = parent(W)
+    (w.a && S.a && w.b && S.b) || throw(ArgumentError())
+    D = _BandedMatrix((-2*(1:∞))', ∞, 1,-1)
+    Mul(Legendre(), D)
+end

src/bases/splines.jl

@@ -58,14 +58,6 @@ function copyto!(dest::SymTridiagonal, AB::Mul2{<:Any,<:Any,<:QuasiAdjoint{<:Any
 end
 
 ## Sub-bases
-# we represent as a Mul with a banded matrix
-function materialize(V::SubQuasiArray{<:Any,2,<:Spline,<:Tuple{<:QSlice,<:AbstractUnitRange}})
-    A = parent(V)
-    _,jr = parentindices(V)
-    first(jr) ≥ 1 || throw(BoundsError())
-    P = BandedMatrix((1-first(jr) => Ones{Int}(length(jr)),), (size(A,2), length(jr)))
-    A*P
-end
 
 
 ## Mass matrix

test/runtests.jl

@@ -1,4 +1,5 @@
-using ContinuumArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, Test
+using ContinuumArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, Test,
+    InfiniteArrays
     import ContinuumArrays: ℵ₁, materialize
     import ContinuumArrays.QuasiArrays: SubQuasiArray
 
@@ -151,3 +152,51 @@ end
 
     @test u[0.1] ≈ 0.00012678835289369413
 end
+
+
+S = Jacobi(true,true)
+W = Diagonal(JacobiWeight(true,true))
+D = Derivative(axes(W,1))
+D*S
+
+D*S[:,1:10]
+
+materialize(D*W*S)
+
+J = Jacobi(false,false)
+J'J
+
+
+S'W'W*S
+
+N = 10
+L = D*W*S[:,1:N]
+    # temporary work around to force 3-term materialize
+    L = *(L.factors[1:3]...) * L.factors[4]
+
+Δ = L'L # weak second derivative
+
+f = Legendre() * Vcat(1:10, Zeros(∞))
+
+
+L'f
+
+axes(L')
+
+Legendre()'f
+
+A = Diagonal(1:∞)
+x = Vcat(1:10, Zeros(∞))
+
+@which A*x
+
+B = BandedMatrices._BandedMatrix((1:∞)', ∞, -1,1)
+
+B*x
+
+N = 100000; @time B[1:N,1:N]
+
+
+
+import LazyArrays: MemoryLayout
+MemoryLayout.(x.arrays)