Revert "Merge branch 'master' into master"

This reverts commit 0bb78e2.

Author: TSGut

src/SemiclassicalOrthogonalPolynomials.jl

@@ -7,14 +7,13 @@ import Base: getindex, axes, size, \, /, *, +, -, summary, ==, copy, sum, unsafe
 
 import ArrayLayouts: MemoryLayout, ldiv, diagonaldata, subdiagonaldata, supdiagonaldata
 import BandedMatrices: bandwidths, AbstractBandedMatrix, BandedLayout, _BandedMatrix
-import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, CachedAbstractVector, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector, ApplyArray,
+import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, CachedAbstractVector, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector,
                     AccumulateAbstractVector, LazyVector, AbstractCachedMatrix, BroadcastLayout
 import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight, NormalizedOPLayout,
                                         Bidiagonal, Tridiagonal, SymTridiagonal, symtridiagonalize, normalizationconstant, LanczosPolynomial,
                                         OrthogonalPolynomialRatio, Weighted, WeightLayout, UnionDomain, oneto, WeightedBasis, HalfWeighted,
                                         golubwelsch, AbstractOPLayout, weight, WeightedOPLayout
 import SingularIntegrals: Associated, associated, Hilbert
-
 import InfiniteArrays: OneToInf, InfUnitRange
 import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayout, MappedBasisLayout, grid, plotgrid, _equals, ExpansionLayout
 import FillArrays: SquareEye
@@ -95,13 +94,20 @@ RaisedOP(Q, y::Number) = RaisedOP(Q, OrthogonalPolynomialRatio(Q,y))
 function jacobimatrix(P::RaisedOP{T}) where T
     ℓ = P.ℓ
     X = jacobimatrix(P.Q)
-    a,b = view(X,band(0)), view(X,band(1)) # a,b = diagonaldata(X), supdiagonaldata(X)
+    a,b = diagonaldata(X), supdiagonaldata(X)
     # non-normalized lower diag of Jacobi
     v = Vcat(zero(T),b .* ℓ)
     c = BroadcastVector((ℓ,a,b,sa,v) -> ℓ*a + b - b*ℓ^2 - sa*ℓ + ℓ*v, ℓ, a, b, a[2:∞], v)
     Tridiagonal(c, BroadcastVector((ℓ,a,b,v) -> a - b * ℓ + v, ℓ,a,b,v), b)
 end
 
+
+
+
+
+
+
+
 """
   the bands of the Jacobi matrix
 """
@@ -136,6 +142,7 @@ SemiclassicalJacobiBand{dv}(a,b,ℓ) where dv = SemiclassicalJacobiBand{dv,promo
 copy(r::SemiclassicalJacobiBand) = r # immutable
 
 
+
 """
    SemiclassicalJacobi(t, a, b, c)
 
@@ -170,7 +177,7 @@ function semiclassical_jacobimatrix(t, a, b, c)
     if a == 0 && b == 0 && c == -1
         # for this special case we can generate the Jacobi operator explicitly
         N = (1:∞)
-        α = neg1c_αcfs(t) # cached coefficients
+        α = T.(neg1c_αcfs(t)) # cached coefficients
         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]./α)
@@ -180,15 +187,10 @@ function semiclassical_jacobimatrix(t, a, b, c)
     P = jacobi(b, a, UnitInterval{T}())
     iszero(c) && return symtridiagonalize(jacobimatrix(P))
     x = axes(P,1)
-    if iseven(c) # QR case
-        return qr_jacobimatrix(x->(t.-x).^(c÷2), Normalized(jacobi(a,b,UnitInterval{T}())))
-    elseif isodd(c)
-        return semiclassical_jacobimatrix(SemiclassicalJacobi(t, a, b, c-1), a, b, c)
-    else # if c is not an even integer, use Lanczos for now
-        return jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), P))
-    end
+    jacobimatrix(LanczosPolynomial(@.(x^a * (1-x)^b * (t-x)^c), P))
 end
 
+
 function symraised_jacobimatrix(Q, y)
     ℓ = OrthogonalPolynomialRatio(Q,y)
     X = jacobimatrix(Q)
@@ -197,35 +199,9 @@ function symraised_jacobimatrix(Q, y)
 end
 
 function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
-    # special cases 
-    if iszero(a) && iszero(b) && c == -1 # (a,b,c) = (0,0,-1) special case
+    if  iszero(a) && iszero(b) && c == -1 # (a,b,c) = (0,0,-1) special case
         semiclassical_jacobimatrix(Q.t, a, b, c)
-    elseif iszero(c) # classical Jacobi polynomial special case
-        jacobimatrix(SemiclassicalJacobi(Q.t, Q.a, Q.b, Q.c))
-    elseif a == Q.a && b == Q.b && c == Q.c # same basis
-        jacobimatrix(Q)
-    end
-
-    # QR decomposition approach for raising
-    Δa = a-Q.a
-    Δb = b-Q.b
-    Δc = c-Q.c
-    if iseven(a-Q.a) && iseven(b-Q.b) && iseven(c-Q.c) && Δa≥0 && Δb≥0 && Δc≥0
-        M = Diagonal(Ones(∞))
-        if !iszero(Δa)
-            M = ApplyArray(*,M,(Q.X)^((a-Q.a)÷2))
-        end
-        if !iszero(Δb)
-            M = ApplyArray(*,M,(I-Q.X)^((b-Q.b)÷2))
-        end
-        if !iszero(Δc)
-            M = ApplyArray(*,M,(Q.t*I-Q.X)^((c-Q.c)÷2))
-        end
-        return qr_jacobimatrix(M,Q,false)
-    end
-
-    # Fallback from decomposition methods
-    if a == Q.a+1 && b == Q.b && c == Q.c  # raising by 1
+    elseif a == Q.a+1 && b == Q.b && c == Q.c  # raising by 1
         symraised_jacobimatrix(Q, 0)
     elseif a == Q.a && b == Q.b+1 && c == Q.c
         symraised_jacobimatrix(Q, 1)
@@ -237,18 +213,20 @@ function semiclassical_jacobimatrix(Q::SemiclassicalJacobi, a, b, c)
         symlowered_jacobimatrix(Q, :a)
     elseif a == Q.a && b == Q.b-1 && c == Q.c
         symlowered_jacobimatrix(Q, :b)
-    elseif a > Q.a  # iterative raising by 1
+    elseif a > Q.a  # iterative raising
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a+1, Q.b, Q.c, Q), a, b, c)
     elseif b > Q.b
         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)
-    elseif b < Q.b  # iterative lowering by 1
+    elseif b < Q.b  # iterative lowering
         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)
     elseif a < Q.a 
         semiclassical_jacobimatrix(SemiclassicalJacobi(Q.t, Q.a-1, Q.b, Q.c, Q), a, b, c)
+    elseif a == Q.a && b == Q.b && c == Q.c # same basis
+        jacobimatrix(Q)
     else
         error("Not Implemented")
     end
@@ -314,46 +292,6 @@ function semijacobi_ldiv(P::SemiclassicalJacobi, Q)
     (P \ R) * _p0(R̃) * (R̃ \ Q)
 end
 
-# returns conversion operator from SemiclassicalJacobi P to SemiclassicalJacobi Q.
-function semijacobi_ldiv(Q::SemiclassicalJacobi,P::SemiclassicalJacobi)
-    @assert Q.t ≈ P.t
-    # For polynomial modifications, use Cholesky. Use Lanzos otherwise.
-    M = Diagonal(Ones(∞))
-    (Q == P) && return M
-    Δa = Q.a-P.a
-    Δb = Q.b-P.b
-    Δc = Q.c-P.c
-    if iseven(Δa) && iseven(Δb) && iseven(Δc)
-        if !iszero(Q.a-P.a)
-            M = ApplyArray(*,M,(P.X)^((Q.a-P.a)÷2))
-        end
-        if !iszero(Q.b-P.b)
-            M = ApplyArray(*,M,(I-P.X)^((Q.b-P.b)÷2))
-        end
-        if !iszero(Q.c-P.c)
-            M = ApplyArray(*,M,(Q.t*I-P.X)^((Q.c-P.c)÷2))
-        end
-        K = qr(M).R
-        return ApplyArray(*, Diagonal(sign.(view(K,band(0))).*Fill(abs.(1/K[1]),∞)), K) # match our normalization choice P_0(x) = 1
-    elseif isinteger(Δa) && isinteger(Δb) && isinteger(Δc)
-        if !iszero(Q.a-P.a)
-            M = ApplyArray(*,M,(P.X)^(Q.a-P.a))
-        end
-        if !iszero(Q.b-P.b)
-            M = ApplyArray(*,M,(I-P.X)^(Q.b-P.b))
-        end
-        if !iszero(Q.c-P.c)
-            M = ApplyArray(*,M,(Q.t*I-P.X)^(Q.c-P.c))
-        end
-        K = cholesky(Symmetric(M)).U
-        return ApplyArray(*, Diagonal(Fill(1/K[1],∞)), K) # match our normalization choice P_0(x) = 1
-    else # fallback option for non-integer weight modification
-        R = SemiclassicalJacobi(P.t, mod(P.a,-1), mod(P.b,-1), mod(P.c,-1))
-        R̃ = toclassical(R)
-        return (P \ R) * _p0(R̃) * (R̃ \ Q)
-    end
-end
-
 struct SemiclassicalJacobiLayout <: AbstractOPLayout end
 MemoryLayout(::Type{<:SemiclassicalJacobi}) = SemiclassicalJacobiLayout()
 
@@ -385,7 +323,6 @@ function copy(L::Ldiv{SemiclassicalJacobiLayout,SemiclassicalJacobiLayout})
     (inv(M_Q) * L') * M_P
 end
 
-\(A::SemiclassicalJacobi, B::SemiclassicalJacobi) = semijacobi_ldiv(A, B)
 \(A::LanczosPolynomial, B::SemiclassicalJacobi) = semijacobi_ldiv(A, B)
 \(A::SemiclassicalJacobi, B::LanczosPolynomial) = semijacobi_ldiv(A, B)
 function \(w_A::WeightedSemiclassicalJacobi, w_B::WeightedSemiclassicalJacobi)
@@ -519,27 +456,7 @@ function SemiclassicalJacobiFamily{T}(data::Vector, t, a, b, c) where T
 end
 
 SemiclassicalJacobiFamily(t, a, b, c) = SemiclassicalJacobiFamily{float(promote_type(typeof(t),eltype(a),eltype(b),eltype(c)))}(t, a, b, c)
-function SemiclassicalJacobiFamily{T}(t, a, b, c) where T
-    ℓa = length(a)
-    ℓb = length(b)
-    ℓc = length(c)
-    # We need to start with a hierarchy containing two entries
-    if  ℓa > 1 && ℓb > 1 && ℓc > 1 
-        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, a[2], b[2], c[2])], t, a, b, c)
-    elseif ℓb > 1 && ℓc > 1 
-        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, first(a), b[2], c[2])], t, a, b, c)
-    elseif ℓa > 1 && ℓc > 1 
-        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, a[2], first(b), c[2])], t, a, b, c)
-    elseif ℓa > 1 && ℓb > 1 
-        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, a[2], b[2], first(c))], t, a, b, c)
-    elseif ℓc > 1 
-        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, first(a), first(b), c[2])], t, a, b, c)
-    elseif ℓa > 1 
-        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, a[2], first(b), first(c))], t, a, b, c)
-    else #if ℓb > 1 
-        return SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c)),SemiclassicalJacobi{T}(t, first(a), b[2], first(c))], t, a, b, c)
-    end
-end
+SemiclassicalJacobiFamily{T}(t, a, b, c) where T = SemiclassicalJacobiFamily{T}([SemiclassicalJacobi{T}(t, first(a), first(b), first(c))], t, a, b, c)
 
 Base.broadcasted(::Type{SemiclassicalJacobi}, t::Number, a::Number, b::Number, c::Number) = SemiclassicalJacobi(t, a, b, c)
 Base.broadcasted(::Type{SemiclassicalJacobi{T}}, t::Number, a::Number, b::Number, c::Number) where T = SemiclassicalJacobi{T}(t, a, b, c)
@@ -555,7 +472,7 @@ _broadcast_getindex(a::Number,k) = a
 function LazyArrays.cache_filldata!(P::SemiclassicalJacobiFamily, inds::AbstractUnitRange)
     t,a,b,c = P.t,P.a,P.b,P.c
     for k in inds
-        P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), P.data[k-2])
+        P.data[k] = SemiclassicalJacobi(t, _broadcast_getindex(a,k), _broadcast_getindex(b,k), _broadcast_getindex(c,k), P.data[k-1])
     end
     P
 end

src/lowering.jl

@@ -1,5 +1,5 @@
 ###########
-# Generic fallback methods for obtaining Jacobi matrices with one of the a, b and c parameters lowered by 1.
+# Generic methods for obtaining Jacobi matrices with one of the a, b and c parameters lowered by 1.
 # passing symbols :a, :b or :c into lowindex determines which parameter is lowered
 ######
 function initialα_gen(P::SemiclassicalJacobi, lowindex::Symbol)
@@ -163,8 +163,8 @@ end
 
 ###########
 # Methods for the special case of computing SemiclassicalJacobi(t,0,0,-1)
-# As this can be built from SemiclassicalJacobi(t,0,0,0) which is just shifted and normalized Legendre(), we have more explicit methods at our disposal due to explicitly known coefficients for the Legendre bases.
 ######
+# As this can be built from SemiclassicalJacobi(t,0,0,0) which is just shifted and normalized Legendre(), we have more explicit methods at our disposal due to explicitly known coefficients for the Legendre bases.
 
 ###
 #   α coefficients implementation, cached, direct and via recurrences

src/twoband.jl

@@ -98,6 +98,7 @@ function jacobimatrix(R::TwoBandJacobi{T}) where T
     # M = (L / L[1,1])' # equal to R.Q \ R.P
 
     Tridiagonal(Interlace(L.dv/L[1,1], (ρ^2-1) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv,L.ev/(-L[1,1])))
+    # Tridiagonal(Interlace(L.dv/L[1,1], (1-ρ^2) * L.ev), Zeros{T}(∞), Interlace((1-ρ^2) * L.dv, L.ev/L[1,1]))
 end
 
 

test/runtests.jl

@@ -117,7 +117,7 @@ end
         @testset "Mass matrix" begin
             @test (T'*(w_T .* T))[1:10,1:10] ≈ sum(w_T)I
             M = W'*(w_W .* W)
-            @test (sum(w_T)*inv(R[1:10,1:10])'L[1:10,1:10]) ≈ M[1:10,1:10]
+            @test (sum(w_T)*inv(R)'L)[1:10,1:10] ≈ M[1:10,1:10]
             @test T[0.1,1:10]' ≈ W[0.1,1:10]' * R[1:10,1:10]
         end
     end
@@ -450,4 +450,3 @@ end
 include("test_derivative.jl")
 include("test_twoband.jl")
 include("test_lowering.jl")
-