Layout issue for conversion with `b=-1`

[DanielVandH]

SemiclassicalOrthogonalPolynomials.jl/src/SemiclassicalOrthogonalPolynomials.jl

Lines 453 to 463 in ed25acf

	elseif A.b == B.b == -1
	Bᵗᵃ¹ᶜ = SemiclassicalJacobi(B.t, B.a, one(B.b), B.c, B)
	Aᵗᵃ¹ᶜ = SemiclassicalJacobi(A.t, A.a, one(A.b), A.c, A)
	Rₐ₁ᵪᵗᵘ¹ᵛ = Aᵗᵃ¹ᶜ \ Bᵗᵃ¹ᶜ
	# Make 1 ⊕ Rₐ₁ᵪᵗᵘ¹ᵛ
	V = eltype(Rₐ₁ᵪᵗᵘ¹ᵛ)
	Rₐ₋₁ᵪᵗᵘ⁻¹ᵛ = Vcat(
	Hcat(one(V), Zeros{V}(1, ∞)),
	Hcat(Zeros{V}(∞), Rₐ₁ᵪᵗᵘ¹ᵛ)
	)
	return Rₐ₋₁ᵪᵗᵘ⁻¹ᵛ

I'm trying to find a better way to do the above conversion since it's giving me some layout headaches. Basically, how can I better compute this conversion matrix of the form $\textrm{diag}(1, \boldsymbol R)$ with better layouts (e.g. colsupport should give the expected results and it should be a banded matrix if R is) than this awkward Vcat/Hcat mess? I could maybe rewrap it as a banded matrix but it might not necessarily be banded e.g. $\boldsymbol R$ could be a dense upper triangular matrix.

To see the issues that I'm running into, consider

using SemiclassicalOrthogonalPolynomials, FillArrays, ArrayLayouts, BandedMatrices
P = SemiclassicalJacobi(2.0, -1 / 2, -1.0, -1 / 2)
Q = SemiclassicalJacobi(2.0, 1 / 2, -1.0, 1 / 2)
R = Q \ P

julia> R
vcat(hcat(Float64, 1×ℵ₀ Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(1)×OneToInf(), hcat(ℵ₀-element Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}} with indices OneToInf(), (ℵ₀×ℵ₀ UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, Float64}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}, BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}}}, Float64}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
 1.0   ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅        …
  ⋅   1.0  0.879719  -0.16154     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅   1.03167    0.867507  -0.172406    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅    ⋅         1.0391     0.866366  -0.175603    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅    ⋅          ⋅         1.04182    0.866144  -0.176996    ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅    ⋅          ⋅          ⋅         1.0431     0.866077  -0.177733    ⋅          ⋅          ⋅          ⋅        …
  ⋅    ⋅    ⋅          ⋅          ⋅          ⋅         1.0438     0.866052  -0.178171    ⋅          ⋅          ⋅
  ⋅    ⋅    ⋅          ⋅          ⋅          ⋅          ⋅         1.04422    0.86604   -0.178453    ⋅          ⋅
  ⋅    ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04449    0.866034  -0.178645    ⋅
  ⋅    ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04468    0.866031  -0.178782
  ⋅    ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04481    0.866029  …
  ⋅    ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04491
  ⋅    ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
 ⋮                                          ⋮                                                      ⋮                    ⋱

Then

julia> colsupport(R, 10) # could be 8:10
1:10

julia> colsupport(R, 5) # could be 3:5
1:5

This seems to lead to issues when multiplying R by other banded matrices. For example,

D = BandedMatrices._BandedMatrix(Zeros(1, ∞), ∞, 0, 0)
R * D

julia> R * D
(vcat(hcat(Float64, 1×ℵ₀ Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(1)×OneToInf(), hcat(ℵ₀-element Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}} with indices OneToInf(), (ℵ₀×ℵ₀ UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, Float64}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}, BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}}}, Float64}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (0, 0) with data 1×ℵ₀ Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf() with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                        ⋮                        ⋮                        ⋮                        ⋮                        ⋱

julia> colsupport(R*D, 1)
1:∞

and e.g.

julia> (R*D)*D
(vcat(hcat(Float64, 1×ℵ₀ Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf()) with indices Base.OneTo(1)×OneToInf(), hcat(ℵ₀-element Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}} with indices OneToInf(), (ℵ₀×ℵ₀ UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, Float64}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}, BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}}}, Float64}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (0, 0) with data 1×ℵ₀ Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf() with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ BandedMatrix{Float64} with bandwidths (0, 0) with data 1×ℵ₀ Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}} with indices Base.OneTo(1)×OneToInf() with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
Error showing value of type ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{Float64, 2, typeof(vcat), Tuple{ApplyArray{Float64, 2, typeof(hcat), Tuple{Float64, Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}}}}, ApplyArray{Float64, 2, typeof(hcat), Tuple{Zeros{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}, ApplyArray{Float64, 2, typeof(*), Tuple{UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, Float64}}}, UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}, BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}}}, Float64}}}}}}}}}, BandedMatrix{Float64, Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}}, InfiniteArrays.OneToInf{Int64}}, BandedMatrix{Float64, Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, InfiniteArrays.OneToInf{Int64}}}, InfiniteArrays.OneToInf{Int64}}}}:
ERROR: MethodError: no method matching materialize!(::MulAdd{ScalarLayout, BandedMatrices.BandedColumns{…}, DualLayout{…}, Float64, Float64, SubArray{…}, Transpose{…}})
The function `materialize!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  materialize!(::Any, ::Base.Broadcast.Broadcasted{BS}) where {NDims, BS<:BlockArrays.AbstractBlockStyle{NDims}}
   @ BlockArrays C:\Users\djv23\.julia\packages\BlockArrays\VccC2\src\blockbroadcast.jl:126
  materialize!(::Any, ::Base.Broadcast.Broadcasted{<:Any})
   @ Base broadcast.jl:874
  materialize!(::Any, ::Any)
   @ Base broadcast.jl:870
  ...

Stacktrace:
  [1] muladd!(α::Float64, A::Float64, B::SubArray{Float64, 2, BandedMatrix{Float64, Zeros{Float64, 2, Tuple{…}}, Base.OneTo{Int64}}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, false}, β::Float64, C::Transpose{Float64, Vector{Float64}}; kw::@Kwargs{})    
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\48qDX\src\muladd.jl:75
  [2] materialize!(M::MulAdd{DualLayout{…}, BandedMatrices.BandedColumns{…}, DualLayout{…}, Float64, Transpose{…}, BandedMatrix{…}, Transpose{…}})
    @ LazyArraysBandedMatricesExt C:\Users\djv23\.julia\packages\LazyArrays\zl2b5\ext\LazyArraysBandedMatricesExt.jl:198
  [3] muladd!(α::Float64, A::Transpose{Float64, ApplyArray{…}}, B::BandedMatrix{Float64, Zeros{…}, Base.OneTo{…}}, β::Float64, C::Transpose{Float64, ApplyArray{…}}; kw::@Kwargs{Czero::Bool})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\48qDX\src\muladd.jl:75
  [4] copyto!
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\48qDX\src\muladd.jl:82 [inlined]
  [5] copy(M::MulAdd{DualLayout{LazyArrays.PaddedRows{…}}, BandedMatrices.BandedColumns{ZerosLayout}, DualLayout{ZerosLayout}, Float64, Transpose{Float64, ApplyArray{…}}, BandedMatrix{Float64, Zeros{…}, Base.OneTo{…}}, Transpose{Float64, Zeros{…}}})       
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\48qDX\src\muladd.jl:77
  [6] copy(M::Mul{DualLayout{LazyArrays.PaddedRows{UnknownLayout}}, BandedMatrices.BandedColumns{ZerosLayout}, Transpose{Float64, ApplyArray{Float64, 1, typeof(vcat), Tuple{…}}}, BandedMatrix{Float64, Zeros{Float64, 2, Tuple{…}}, Base.OneTo{Int64}}})      
    @ LazyArraysBandedMatricesExt C:\Users\djv23\.julia\packages\LazyArrays\zl2b5\ext\LazyArraysBandedMatricesExt.jl:635
  [7] materialize
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\48qDX\src\mul.jl:137 [inlined]
  [8] mul
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\48qDX\src\mul.jl:138 [inlined]
  [9] *(A::Transpose{Float64, ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}}}}, B::BandedMatrix{Float64, Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Base.OneTo{Int64}})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\48qDX\src\mul.jl:189
 [10] *(tu::Transpose{Float64, ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}}}}, B::BandedMatrix{Float64, Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Base.OneTo{Int64}}, v::Vector{Float64})
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\matmul.jl:1115
 [11] _transposefirst_andmul(::ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, Zeros{Float64, 1, Tuple{Base.OneTo{Int64}}}}}, ::BandedMatrix{Float64, Zeros{Float64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Base.OneTo{Int64}}, ::Vararg{Any})    
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\zl2b5\src\linalg\mul.jl:211
 [12] _mul_getindex(args::Tuple{ApplyArray{Float64, 2, typeof(vcat), Tuple{…}}, BandedMatrix{Float64, Zeros{…}, InfiniteArrays.OneToInf{…}}, BandedMatrix{Float64, Zeros{…}, InfiniteArrays.OneToInf{…}}}, k::Int64, j::Int64)
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\zl2b5\src\linalg\mul.jl:218
 [13] getindex
    @ C:\Users\djv23\.julia\packages\LazyArrays\zl2b5\src\linalg\mul.jl:184 [inlined]
 [14] isassigned(::ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{Float64, 2, typeof(vcat), Tuple{…}}, BandedMatrix{Float64, Zeros{…}, InfiniteArrays.OneToInf{…}}, BandedMatrix{Float64, Zeros{…}, InfiniteArrays.OneToInf{…}}}}, ::Int64, ::Int64)        
    @ Base .\multidimensional.jl:1612
 [15] alignment(io::IOContext{Base.TTY}, X::AbstractVecOrMat, rows::Vector{Int64}, cols::Vector{Int64}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64, ncols::Infinities.InfiniteCardinal{0})
    @ Base .\arrayshow.jl:68
 [16] _print_matrix(io::IOContext{Base.TTY}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::InfiniteArrays.InfUnitRange{Int64}, colsA::InfiniteArrays.InfUnitRange{Int64})
    @ Base .\arrayshow.jl:207
 [17] print_matrix(io::IOContext{Base.TTY}, X::ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{…}, BandedMatrix{…}, BandedMatrix{…}}}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
    @ Base .\arrayshow.jl:171
 [18] print_matrix
    @ .\arrayshow.jl:171 [inlined]
 [19] print_array
    @ .\arrayshow.jl:358 [inlined]
 [20] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, X::ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{…}, BandedMatrix{…}, BandedMatrix{…}}})
    @ Base .\arrayshow.jl:399
 [21] (::REPL.var"#68#69"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:367
 [22] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:661
 [23] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:353
 [24] display(d::REPL.REPLDisplay, x::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:372
 [25] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:340
 [26] #invokelatest#2
    @ .\essentials.jl:1055 [inlined]
 [27] invokelatest
    @ .\essentials.jl:1052 [inlined]
 [28] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:409
 [29] (::REPL.var"#70#71"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:378
 [30] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:661
 [31] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:376
 [32] (::REPL.var"#do_respond#96"{Bool, Bool, REPL.var"#112#130"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:1003
 [33] (::VSCodeServer.var"#103#106"{REPL.var"#do_respond#96"{Bool, Bool, REPL.var"#112#130"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt}})(mi::REPL.LineEdit.MIState, buf::IOBuffer, ok::Bool)
    @ VSCodeServer c:\Users\djv23\.vscode\extensions\julialang.language-julia-1.127.2\scripts\packages\VSCodeServer\src\repl.jl:122
 [34] #invokelatest#2
    @ .\essentials.jl:1055 [inlined]
 [35] invokelatest
    @ .\essentials.jl:1052 [inlined]
 [36] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\LineEdit.jl:2755
 [37] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:1471
 [38] (::REPL.var"#75#81"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL C:\Users\djv23\.julia\juliaup\julia-1.11.1+0.x64.w64.mingw32\share\julia\stdlib\v1.11\REPL\src\REPL.jl:480
Some type information was truncated. Use `show(err)` to see complete types.

[DanielVandH]

I could maybe rewrap it as a banded matrix but it might not necessarily be banded e.g. R could be a dense upper triangular matrix.

Actually this for some reason I'm also not sure about for a temporary workaround when R is banded since the stored matrix for b=1 isn't a banded matrix where I have access to the banded storage:

julia> P = SemiclassicalJacobi(2.0, -1 / 2, 1.0, -1 / 2)
SemiclassicalJacobi with weight x^-0.5 * (1-x)^1.0 * (2.0-x)^-0.5 on 0..1

julia> Q = SemiclassicalJacobi(2.0, 1 / 2, 1.0, 1 / 2)
SemiclassicalJacobi with weight x^0.5 * (1-x)^1.0 * (2.0-x)^0.5 on 0..1

julia> R = Q \ P
(ℵ₀×ℵ₀ UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, Float64}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{Float64, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}, BroadcastVector{Float64, typeof(-), Tuple{ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}}}, Float64}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
 1.0  0.879719  -0.16154     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅        …  
  ⋅   1.03167    0.867507  -0.172406    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅         1.0391     0.866366  -0.175603    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅          ⋅         1.04182    0.866144  -0.176996    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅          ⋅          ⋅         1.0431     0.866077  -0.177733    ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅          ⋅          ⋅          ⋅         1.0438     0.866052  -0.178171    ⋅          ⋅          ⋅          ⋅        …
  ⋅    ⋅          ⋅          ⋅          ⋅          ⋅         1.04422    0.86604   -0.178453    ⋅          ⋅          ⋅
  ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04449    0.866034  -0.178645    ⋅          ⋅
  ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04468    0.866031  -0.178782    ⋅
  ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04481    0.866029  -0.178883
  ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04491    0.866028  …
  ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅         1.04499
  ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
  ⋅    ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅
 ⋮                                                ⋮                                                      ⋮                    ⋱

julia> R.
args
f

[TSGut]

Try

R[band(0)]
R[band(1)]
R[band(2)]

[DanielVandH]

Oh of course, that'll work for that second case thanks.