From 768dab9d51d0c2924533cc9cf445ecfec1f64934 Mon Sep 17 00:00:00 2001 From: Timon Salar Gutleb Date: Tue, 11 Jul 2023 13:57:23 +0100 Subject: [PATCH 1/2] Add more tests for v4.0 --- test/runtests.jl | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1a98745..2ce632a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ using SemiclassicalOrthogonalPolynomials using ClassicalOrthogonalPolynomials, ContinuumArrays, BandedMatrices, QuasiArrays, Test, LazyArrays, FillArrays, LinearAlgebra import BandedMatrices: _BandedMatrix -import SemiclassicalOrthogonalPolynomials: op_lowering, RaisedOP, jacobiexpansion +import SemiclassicalOrthogonalPolynomials: op_lowering, RaisedOP, jacobiexpansion, semijacobi_ldiv_direct import ClassicalOrthogonalPolynomials: recurrencecoefficients, orthogonalityweight, symtridiagonalize @testset "Jacobi" begin @@ -65,6 +65,21 @@ end @test SemiclassicalJacobi(1.013,2,2,2.3)[0.7, 4] ≈ SemiclassicalJacobi{Float64}(1.013,2,2,2.3)[0.7, 4] end +@testset "Test ldiv versus not recommended direct ldiv" begin + # set 1 + P = SemiclassicalJacobi(1.1,0,0,4) + Q = SemiclassicalJacobi(1.1,1,2,7) + R = Q \ P + Ralt = semijacobi_ldiv_direct(Q,P) + @test R[1:20,1:20] ≈ Ralt[1:20,1:20] + # set 2 + P = SemiclassicalJacobi(1.23,4,1,2) + Q = SemiclassicalJacobi(1.23,7,4,6) + R = Q \ P + Ralt = semijacobi_ldiv_direct(Q,P) + @test R[1:20,1:20] ≈ Ralt[1:20,1:20] +end + @testset "Half-range Chebyshev" begin @testset "T and W" begin T = SemiclassicalJacobi(2, -1/2, 0, -1/2) From 9b40d1294c496bbec71c39a7bfe87770a76413dd Mon Sep 17 00:00:00 2001 From: Timon Salar Gutleb Date: Tue, 11 Jul 2023 14:09:24 +0100 Subject: [PATCH 2/2] tidy up comments somewhat --- src/SemiclassicalOrthogonalPolynomials.jl | 22 +++++++++------------- src/neg1c.jl | 6 ++---- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/src/SemiclassicalOrthogonalPolynomials.jl b/src/SemiclassicalOrthogonalPolynomials.jl index 7fec388..0dc1580 100644 --- a/src/SemiclassicalOrthogonalPolynomials.jl +++ b/src/SemiclassicalOrthogonalPolynomials.jl @@ -128,19 +128,16 @@ SemiclassicalJacobi{T}(t, a, b, c) where T = SemiclassicalJacobi(convert(T,t), c WeightedSemiclassicalJacobi{T}(t, a, b, c, P...) where T = SemiclassicalJacobiWeight{T}(convert(T,t), convert(T,a), convert(T,b), convert(T,c)) .* SemiclassicalJacobi{T}(convert(T,t), convert(T,a), convert(T,b), convert(T,c), P...) - - function semiclassical_jacobimatrix(t, a, b, c) T = float(promote_type(typeof(t), typeof(a), typeof(b), typeof(c))) if iszero(a) && iszero(b) && c == -one(T) # for this special case we can generate the Jacobi operator explicitly N = (1:∞) - α = neg1c_αcfs(one(T)*t) # cached coefficients + α = neg1c_αcfs(one(T)*t) A = Vcat((α[1]+1)/2 , -N./(N.*4 .- 2).*α .+ (N.+1)./(N.*4 .+ 2).*α[2:end].+1/2) C = -(N)./(N.*4 .- 2) B = Vcat((α[1]^2*3-α[1]*α[2]*2-1)/6 , -(N)./(N.*4 .+ 2).*α[2:end]./α) - # if J is Tridiagonal(c,a,b) then for norm. OPs it becomes SymTridiagonal(a, sqrt.( b.* c)) - return SymTridiagonal(A, sqrt.(B.*C)) + return SymTridiagonal(A, sqrt.(B.*C)) # if J is Tridiagonal(c,a,b) then for norm. OPs it becomes SymTridiagonal(a, sqrt.( b.* c)) end P = Normalized(jacobi(b, a, UnitInterval{T}())) iszero(c) && return jacobimatrix(P) @@ -150,7 +147,7 @@ function semiclassical_jacobimatrix(t, a, b, c) return qr_jacobimatrix(x->(t-x),P) elseif isinteger(c) && c ≥ 0 # reduce other integer c cases to hierarchy return SemiclassicalJacobi.(t, a, b, 0:Int(c))[end].X - else # if c is not an integer, use Lanczos for now + else # if c is not an integer, use Lanczos x = axes(P,1) return jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), jacobi(b, a, UnitInterval{T}()))) end @@ -188,7 +185,7 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c) semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b+1, Q.c, Q), a, b, c) elseif c > Q.c semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c+1, Q), a, b, c) - # TODO: Implement lowering via QL/reverse Cholesky or via inverting R + # TODO: Implement lowering via QL/reverse Cholesky or via inverting R? # elseif b < Q.b # iterative lowering by 1 # semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b-1, Q.c, Q), a, b, c) # elseif c < Q.c @@ -233,8 +230,7 @@ Gives the Lowering operator from OPs w.r.t. (x-y)*w(x) to Q as constructed from Chistoffel–Darboux """ function op_lowering(Q, y) - # we first use Christoff-Darboux with d = 1 - # But we want the first OP to be 1 so we rescale + # we first use Christoff-Darboux with d = 1 but the first OP should be 1 so we rescale P = RaisedOP(Q, y) A,_,_ = recurrencecoefficients(Q) d = -inv(A[1]*_p0(Q)*P.ℓ[1]) @@ -268,16 +264,16 @@ function semijacobi_ldiv_direct(Q::SemiclassicalJacobi, P::SemiclassicalJacobi) M = Diagonal(Ones(∞)) if iseven(Δa) && iseven(Δb) && iseven(Δc) M = qr(P.X^(Δa÷2)*(I-P.X)^(Δb÷2)*(Q.t*I-P.X)^(Δc÷2)).R - return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) # match normalization choice P_0(x) = 1 + return ApplyArray(*, Diagonal(sign.(view(M,band(0))).*Fill(abs.(1/M[1]),∞)), M) elseif isone(Δa) && iszero(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (1,0,0) M = cholesky(P.X).U - return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1 + return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) elseif iszero(Δa) && isone(Δb) && iszero(Δc) # special case (Δa,Δb,Δc) = (0,1,0) M = cholesky(I-P.X).U - return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1 + return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) elseif iszero(Δa) && iszero(Δb) && isone(Δc) # special case (Δa,Δb,Δc) = (0,0,1) M = cholesky(Q.t*I-P.X).U - return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1 + return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc) M = cholesky(Symmetric(P.X^(Δa)*(I-P.X)^(Δb)*(Q.t*I-P.X)^(Δc))).U return ApplyArray(*, Diagonal(Fill(1/M[1],∞)), M) # match normalization choice P_0(x) = 1 diff --git a/src/neg1c.jl b/src/neg1c.jl index 3cb7d3f..efae617 100644 --- a/src/neg1c.jl +++ b/src/neg1c.jl @@ -41,7 +41,7 @@ function backαcoeff!(α, t, inds) end end -# cached implementation of normalization constants +# cached implementation of normalization constants mutable struct neg1c_normconstant{T} <: AbstractCachedVector{T} data::Vector{T} t::T @@ -79,9 +79,7 @@ function cache_filldata!(B::neg1c_normconstant{T}, inds::UnitRange{Int}) where T B.data[m+1:n] = norm[1:end-1] end -### -# bidiagonal operator which converts from OPs wrt (t-x)^-1 to shifted Legendre -### +# bidiagonal operator which converts from OPs wrt (t-x)^-1 to shifted Legendre function neg1c_tolegendre(t::T) where T nc = neg1c_normconstant(t) α = neg1c_αcfs(t)