big steps towards generic OP basis

Author: TSGut

Project.toml

@@ -22,6 +22,7 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+SemiclassicalOrthogonalPolynomials = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]

src/decompOPs.jl

@@ -1,26 +1,25 @@
 # This currently takes the weight multiplication operator as input.
 # I will probably change this to take the weight function instead.
-function cholesky_jacobimatrix(W::Symmetric)
-    bands = CholeskyJacobiBands(W) # the cached array only needs to store two bands bc of symmetry
+function cholesky_jacobimatrix(w, P::OrthogonalPolynomial)
+    W = Symmetric(P \ (w.(axes(P,1)) .* P)) # Compute weight multiplication via Clenshaw
+    bands = CholeskyJacobiBands(W, P) # the cached array only needs to store two bands bc of symmetry
     return SymTridiagonal(bands[1,:],bands[2,:])
 end
 
 # The generated Jacobi operators are symmetric tridiagonal, so we store their data as two bands
 mutable struct CholeskyJacobiBands{T} <: AbstractCachedMatrix{T}
     data::Matrix{T}
     U::UpperTriangular
-    X::Symmetric{T}
+    X::SymTridiagonal{T}
     datasize::Int
     array
 end
 
-# SymTridiagonal currently doesn't parse as Symmetric, so here's a Q&D workaround for conversion
-symmjacobim(J::SymTridiagonal) = Symmetric(BandedMatrix(0=>J.dv, 1=>J.ev))
-
 # Computes the initial data for the Jacobi operator bands
-function CholeskyJacobiBands(W::Symmetric{T}) where T
+function CholeskyJacobiBands(W::Symmetric{T}, P::OrthogonalPolynomial) where T
+    @assert P isa Normalized
     U = cholesky(W).U
-    X = symmjacobim(jacobimatrix(Normalized(Legendre()[affine(zero(T)..one(T),Inclusion(-one(T)..one(T))),:])))
+    X = jacobimatrix(P)
     dat = zeros(T,2,10)
     for k in 1:10
         dat[1,k] = (U * (X * (U \ [zeros(k-1); 1; zeros(∞)])))[k]