Default recurrencecoefficients (#27)

* test adjoint views work

* recurrencecoefficients now defaults to using jacobimatrix

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.3.3"
+version = "0.3.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -28,7 +28,7 @@ ArrayLayouts = "0.6.2"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.14.1, 0.15"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.6.2"
+ContinuumArrays = "0.6.4"
 DomainSets = "0.4, 0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"

src/ClassicalOrthogonalPolynomials.jl

@@ -144,7 +144,11 @@ The relationship with the Jacobi matrix is:
 C[n+1]/A[n+1] == X[n,n+1]
 ```
 """
-recurrencecoefficients(P) = error("Override for $(typeof(P))")
+function recurrencecoefficients(Q::AbstractQuasiMatrix{T}) where T
+    X = jacobimatrix(Q)
+    c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
+    inv.(c), -(a ./ c), Vcat(zero(T), b) ./ c
+end
 
 
 const WeightedOrthogonalPolynomial{T, A<:AbstractQuasiVector, B<:OrthogonalPolynomial} = WeightedBasis{T, A, B}

src/lanczos.jl

@@ -183,7 +183,6 @@ axes(Q::LanczosPolynomial) = (axes(Q.w,1),OneToInf())
 _p0(Q::LanczosPolynomial) = inv(Q.data.γ[1])*_p0(Q.P)
 
 jacobimatrix(Q::LanczosPolynomial) = SymTridiagonal(LanczosJacobiBand(Q.data, :d), LanczosJacobiBand(Q.data, :du))
-recurrencecoefficients(Q::LanczosPolynomial) = normalized_recurrencecoefficients(Q)
 
 Base.summary(io::IO, Q::LanczosPolynomial{T}) where T = print(io, "LanczosPolynomial{$T} with weight with singularities $(singularities(Q.w))")
 Base.show(io::IO, Q::LanczosPolynomial{T}) where T = summary(io, Q)

src/normalized.jl

@@ -55,17 +55,6 @@ axes(Q::Normalized) = axes(Q.P)
 
 _p0(Q::Normalized) = Q.scaling[1]
 
-# p_{n+1} = (A_n * x + B_n) * p_n - C_n * p_{n-1}
-# q_{n+1}/h[n+1] = (A_n * x + B_n) * q_n/h[n] - C_n * p_{n-1}/h[n-1]
-# q_{n+1} = (h[n+1]/h[n] * A_n * x + h[n+1]/h[n] * B_n) * q_n - h[n+1]/h[n-1] * C_n * p_{n-1}
-
-function normalized_recurrencecoefficients(Q::AbstractQuasiMatrix{T}) where T
-    X = jacobimatrix(Q)
-    a,b = diagonaldata(X), supdiagonaldata(X)
-    inv.(b), -(a ./ b), Vcat(zero(T), b) ./ b
-end
-
-recurrencecoefficients(Q::Normalized{T}) where T = normalized_recurrencecoefficients(Q)
 
 # x * p[n] = c[n-1] * p[n-1] + a[n] * p[n] + b[n] * p[n+1]
 # x * q[n]/h[n] = c[n-1] * q[n-1]/h[n-1] + a[n] * q[n]/h[n] + b[n] * q[n+1]/h[n+1]

test/test_chebyshev.jl

@@ -349,6 +349,13 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
         @test chebyshevu.(0:5, BigFloat(1)/10) == ChebyshevU{BigFloat}()[BigFloat(1)/10, 1:6]
         @test chebyshevu.(0:5, 2) == Base.unsafe_getindex(ChebyshevU(), 2, 1:6)
     end
+
+    @testset "Adjoint views" begin
+        T = ChebyshevT(); U = ChebyshevU()
+        D = Derivative(axes(T,1))
+        V = view(T,:,[1,3,4])
+        @test (U\(D*V))[1:5,:] == (U \ (V'D')')[1:5,:] == (U\(D*T))[1:5,[1,3,4]] 
+    end
 end
 
 struct QuadraticMap{T} <: Map{T} end