search

JuliaApproximation / SemiclassicalOrthogonalPolynomials.jl

A Julia repository for semiclassical orthogonal polynomials
MIT License
7 stars 3 forks source link

Can't invert connection matrices #108

Closed DanielVandH closed 6 days ago

DanielVandH commented 1 month ago 
julia> using SemiclassicalOrthogonalPolynomials, ArrayLayouts

julia> R = SemiclassicalJacobi(2.0, 1.0, 1.0, 1.0) \ SemiclassicalJacobi(2.0, 1.0, 0.0, 1.0);

julia> inv(R)
ERROR: MethodError: no method matching Matrix{…}(::LinearAlgebra.UniformScaling{…}, ::Tuple{…})

Closest candidates are:
  Array{T, N}(::Missing, ::Any...) where {T, N}
   @ Base baseext.jl:43
  Array{T, N}(::Nothing, ::Any...) where {T, N}
   @ Base baseext.jl:42
  Matrix{T}(::UndefInitializer, ::Tuple{Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}, Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}}) where T
   @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:7
  ...

Stacktrace:
 [1] inv(A::LinearAlgebra.UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{…}})
   @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:797
 [2] map
   @ .\tuple.jl:292 [inlined]
 [3] _inv(Lay::LazyArrays.ApplyLayout{…}, ::Tuple{…}, A::LazyArrays.ApplyArray{…})
   @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:392
 [4] inv(A::LazyArrays.ApplyArray{Float64, 2, typeof(*), Tuple{LinearAlgebra.Diagonal{…}, LinearAlgebra.UpperTriangular{…}}})
   @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\factorizations.jl:427
 [5] top-level scope
   @ REPL[19]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> MemoryLayout(R), map(MemoryLayout, R.args)
(LazyArrays.ApplyLayout{typeof(*)}(), (DiagonalLayout{FillLayout}(), TriangularLayout{'U', 'N', InfiniteLinearAlgebra.AdaptiveLayout{BandedMatrices.BandedColumns{DenseColumnMajor}}}()))

I guess R should have something like a memory layout of ApplyBandedLayout{typeof(*)}() instead? I'm currently just doing ApplyArray(inv, R)

TSGut commented 1 month ago

You could also consider ApplyArray(\,R,Eye{T}(∞)) instead of ApplyArray(inv, R). If I remember correctly I had to use this somewhere and the former behaved better.

DanielVandH commented 1 month ago

Hm. Using that approach gives me errors down the line related to trying to create infinite matrices (meaning ArgumentError: Cannot create infinite Array) for what I'm working with. Do you know where you used it previously?

DanielVandH commented 1 month ago

I am able to get the ApplyBandedLayout that I mentioned using 

using ClassicalOrthogonalPolynomials, SemiclassicalOrthogonalPolynomials, ArrayLayouts
A = SemiclassicalJacobi(2.0, 2.0, 1.0, 1.0) \ SemiclassicalJacobi(2.0, 1.0, 0.0, 1.0)
AA = Jacobi(2.0, 1.0) \ Jacobi(1.0, 0.0)

MemoryLayout(A)  # LazyArrays.ApplyLayout{typeof(*)}()
MemoryLayout(AA) # LazyArraysBandedMatricesExt.ApplyBandedLayout{typeof(*)}()

using SemiclassicalOrthogonalPolynomials, ArrayLayouts, InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, LazyArrays, InfiniteArrays
ext = Base.get_extension(LazyArrays, :LazyArraysBandedMatricesExt)
ArrayLayouts.MemoryLayout(::Type{<:UpperTriangular{V, <:InfiniteLinearAlgebra.AdaptiveCholeskyFactors{T, <:BandedMatrices.AbstractBandedMatrix}}}) where {V, T} = BandedMatrices.BandedLayout
LazyArrays.applylayout(::Type{typeof(*)}, ::DiagonalLayout, ::Type{<:ext.BandedLayouts}) = ext.ApplyBandedLayout{typeof(*)}()

MemoryLayout(A) # LazyArraysBandedMatricesExt.ApplyBandedLayout{typeof(*)}()

although the new MemoryLayout method seems a bit clunky as if I should be able to do something better... it's still not enough anyway:

julia> A[1, 1]
ERROR: MethodError: sublayout(::LazyArraysBandedMatricesExt.ApplyBandedLayout{typeof(*)}, ::Type{Tuple{Int64, InfiniteArrays.OneToInf{Int64}}}) is ambiguous.

Candidates:
  sublayout(lay::LazyArraysBandedMatricesExt.ApplyBandedLayout, ind)
    @ LazyArraysBandedMatricesExt C:\Users\User\.julia\packages\LazyArrays\MLFsy\ext\LazyArraysBandedMatricesExt.jl:226
  sublayout(::AbstractBandedLayout, ::Type{<:Tuple{Integer, JR}}) where JR<:Union{Base.IdentityUnitRange{<:InfiniteArrays.AbstractInfUnitRange{<:Integer}}, Base.Slice{<:InfiniteArrays.AbstractInfUnitRange{<:Integer}}, InfiniteArrays.AbstractInfUnitRange{<:Integer}, InfiniteArrays.InfStepRange{<:Integer}}  
    @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:305

Possible fix, define
  sublayout(::LazyArraysBandedMatricesExt.ApplyBandedLayout, ::Type{…}) where JR<:Union{…}

Stacktrace:
 [1] MemoryLayout
   @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\memorylayout.jl:186 [inlined]
 [2] MemoryLayout
   @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\memorylayout.jl:169 [inlined]
 [3] sub_materialize
   @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
 [4] layout_getindex
   @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
 [5] getindex
   @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:159 [inlined]
 [6] map
   @ .\tuple.jl:342 [inlined]
 [7] _mul_getindex(args::Tuple{ApplyArray{…}, ApplyArray{…}}, k::Int64, j::Int64)
   @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:218
 [8] getindex(::ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{…}, ApplyArray{…}}}, ::Int64, ::Int64)
   @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:184
 [9] top-level scope
   @ REPL[5]:1
Some type information was truncated. Use `show(err)` to see complete types.

I'm not sure how to get sublayout defined perfectly. The two candidates give

II = Tuple{Int64, InfiniteArrays.OneToInf{Int64}}
invoke(ArrayLayouts.sublayout, Tuple{AbstractBandedLayout, Any}, MemoryLayout(A), II)
invoke(ArrayLayouts.sublayout, Tuple{ext.ApplyBandedLayout, typeof(II)}, MemoryLayout(A), II)
julia> invoke(ArrayLayouts.sublayout, Tuple{AbstractBandedLayout, Any}, MemoryLayout(A), II)
UnknownLayout()

julia> invoke(ArrayLayouts.sublayout, Tuple{ext.ApplyBandedLayout, typeof(II)}, MemoryLayout(A), II)
LazyArrays.ApplyLayout{typeof(*)}()

so I tried 

ArrayLayouts.sublayout(::ext.ApplyBandedLayout{F}, ::Type{<:Tuple{Integer, JR}}) where {F, JR<:InfiniteArrays.InfAxes} = LazyArrays.ApplyLayout{F}()

but A[1, 1] just hangs.

julia> A[1, 1]
ERROR: InterruptException:
Stacktrace:
  [1] partialcholesky!(F::InfiniteLinearAlgebra.AdaptiveCholeskyFactors{…}, n::Int64)
    @ InfiniteLinearAlgebra C:\Users\User\.julia\packages\InfiniteLinearAlgebra\cI65l\src\infcholesky.jl:25
  [2] getindex
    @ C:\Users\User\.julia\packages\InfiniteLinearAlgebra\cI65l\src\infcholesky.jl:45 [inlined]
  [3] getindex
    @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:241 [inlined]
  [4] getindex
    @ .\subarray.jl:290 [inlined]
  [5] _default_blasmul_loop!
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:164 [inlined]
  [6] default_blasmul!(α::Float64, A::SubArray{…}, B::BandedMatrix{…}, β::Bool, C::Matrix{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:188
  [7] materialize!(M::MulAdd{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:264
  [8] muladd!
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:75 [inlined]
  [9] copyto!
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:82 [inlined]
 [10] copy(M::MulAdd{…})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\muladd.jl:77
 [11] copy
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\mul.jl:140 [inlined]
 [12] materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\mul.jl:137 [inlined]
 [13] mul
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\mul.jl:138 [inlined]
 [14] *(A::SubArray{Float64, 2, UpperTriangular{…}, Tuple{…}, false}, B::BandedMatrix{Float64, Matrix{…}, Base.OneTo{…}})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\mul.jl:227
 [15] _fillcholeskybanddata!(J::ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, inds::UnitRange{Int64})
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\DiRvH\src\choleskyQR.jl:119
 [16] resizedata!(K::ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, n::Int64, m::Int64)
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\DiRvH\src\choleskyQR.jl:107
 [17] getindex(K::ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, n::Int64, m::Int64)
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\DiRvH\src\choleskyQR.jl:95
 [18] getindex
    @ .\subarray.jl:290 [inlined]
 [19] _getindex
    @ .\abstractarray.jl:1341 [inlined]
 [20] getindex
    @ .\abstractarray.jl:1291 [inlined]
 [21] _broadcast_getindex
    @ .\broadcast.jl:662 [inlined]
 [22] _getindex
    @ .\broadcast.jl:706 [inlined]
 [23] _getindex
    @ .\broadcast.jl:705 [inlined]
 [24] _broadcast_getindex
    @ .\broadcast.jl:681 [inlined]
 [25] getindex
    @ .\broadcast.jl:636 [inlined]
 [26] getindex
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\lazybroadcasting.jl:93 [inlined]
 [27] getindex(A::LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{…}, BroadcastVector{…}}, i::Int64, j::Int64)
    @ LazyBandedMatrices C:\Users\User\.julia\packages\LazyBandedMatrices\gqWs9\src\tridiag.jl:306
 [28] getindex
    @ .\subarray.jl:290 [inlined]
 [29] (BandedMatrix{…})(A::SubArray{…}, bnds::Tuple{…})
    @ BandedMatrices C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:214
 [30] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:207 [inlined]
 [31] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:219 [inlined]
 [32] _BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:245 [inlined]
 [33] BandedMatrix
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\banded\BandedMatrix.jl:246 [inlined]
 [34] sub_materialize
    @ C:\Users\User\.julia\packages\BandedMatrices\d4g9c\src\generic\indexing.jl:28 [inlined]
 [35] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:131 [inlined]
 [36] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
 [37] layout_getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
 [38] getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:153 [inlined]
 [39] resizedata!(::BandedMatrices.BandedColumns{…}, ::SymTridiagonalLayout{…}, B::LazyArrays.CachedArray{…}, n::Int64, m::Int64)
    @ LazyArraysBandedMatricesExt C:\Users\User\.julia\packages\LazyArrays\MLFsy\ext\LazyArraysBandedMatricesExt.jl:492
 [40] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\cache.jl:218 [inlined]
 [41] partialcholesky!(F::InfiniteLinearAlgebra.AdaptiveCholeskyFactors{…}, n::Int64)
    @ InfiniteLinearAlgebra C:\Users\User\.julia\packages\InfiniteLinearAlgebra\cI65l\src\infcholesky.jl:28
 [42] _fillcholeskybanddata!(J::ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, inds::UnitRange{Int64})
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\DiRvH\src\choleskyQR.jl:116
 [43] resizedata!(K::ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, n::Int64, m::Int64)
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\DiRvH\src\choleskyQR.jl:107
 [44] getindex(C::LazyBandedMatrices.SymTridiagonal{…}, kr::UnitRange{…}, jr::UnitRange{…})
    @ ClassicalOrthogonalPolynomials C:\Users\User\.julia\packages\ClassicalOrthogonalPolynomials\DiRvH\src\choleskyQR.jl:88
 [45] resizedata!(::BandedMatrices.BandedColumns{…}, ::SymTridiagonalLayout{…}, B::LazyArrays.CachedArray{…}, n::Int64, m::Int64)
    @ LazyArraysBandedMatricesExt C:\Users\User\.julia\packages\LazyArrays\MLFsy\ext\LazyArraysBandedMatricesExt.jl:489
 [46] resizedata!
    @ C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\cache.jl:218 [inlined]
 [47] partialcholesky!(F::InfiniteLinearAlgebra.AdaptiveCholeskyFactors{…}, n::Int64)
    @ InfiniteLinearAlgebra C:\Users\User\.julia\packages\InfiniteLinearAlgebra\cI65l\src\infcholesky.jl:28
 [48] getindex
    @ C:\Users\User\.julia\packages\InfiniteLinearAlgebra\cI65l\src\infcholesky.jl:45 [inlined]
 [49] getindex
    @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\adjtrans.jl:329 [inlined]
 [50] getindex
    @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:237 [inlined]
 [51] getindex
    @ .\subarray.jl:290 [inlined]
 [52] _getindex
    @ .\abstractarray.jl:1336 [inlined]
 [53] getindex
    @ .\abstractarray.jl:1291 [inlined]
 [54] _generic_matvecmul!(C::LazyArrays.CachedArray{…}, tA::Char, A::SubArray{…}, B::Vector{…}, _add::LinearAlgebra.MulAddMul{…})
    @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:745
 [55] generic_matvecmul!
    @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:687 [inlined]
 [56] mul!
    @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:66 [inlined]
 [57] mul!
    @ C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:237 [inlined]
 [58] *(A::SubArray{Float64, 2, LowerTriangular{…}, Tuple{…}, false}, x::Vector{Float64})
    @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:57
 [59] sub_materialize(lay::LazyArrays.ApplyLayout{typeof(*)}, V::SubArray{Float64, 1, ApplyArray{…}, Tuple{…}, false})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:261
 [60] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
 [61] layout_getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
 [62] getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:159 [inlined]
 [63] map
    @ .\tuple.jl:342 [inlined]
 [64] _mul_getindex(args::Tuple{ApplyArray{…}, ApplyArray{…}}, k::Int64, j::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:218
 [65] getindex(::ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{…}, ApplyArray{…}}}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:184
Some type information was truncated. Use `show(err)` to see complete types.

I also tried

ArrayLayouts.sublayout(::ext.ApplyBandedLayout{F}, ::Type{<:Tuple{Integer, JR}}) where {F, JR<:InfiniteArrays.InfAxes} = ext.ApplyBandedLayout{F}()

but

julia> A[1, 1]
ERROR: MethodError: no method matching (Zeros{Float64})(::Infinities.NotANumber)

Closest candidates are:
  (Zeros{T})(::Infinities.Infinity) where T
   @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:46
  (Zeros{T})(::AbstractArray) where T
   @ FillArrays C:\Users\User\.julia\packages\FillArrays\eOEVm\src\FillArrays.jl:320
  (Zeros{T})(::SZ) where {T, N, SZ<:Tuple{Vararg{Any, N}}}
   @ FillArrays C:\Users\User\.julia\packages\FillArrays\eOEVm\src\FillArrays.jl:312
  ...

Stacktrace:
  [1] _padded_sub_materialize(v::SubArray{Float64, 1, ApplyArray{…}, Tuple{…}, false})
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\padded.jl:477
  [2] sub_materialize
    @ C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:396 [inlined]
  [3] sub_materialize(L::LazyArrays.PaddedColumns{UnknownLayout}, V::SubArray{Float64, 1, ApplyArray{…}, Tuple{…}, false})
    @ ArrayLayouts C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:131
  [4] sub_materialize
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:132 [inlined]
  [5] layout_getindex
    @ C:\Users\User\.julia\packages\ArrayLayouts\3byqH\src\ArrayLayouts.jl:138 [inlined]
  [6] _unsafe_getindex(::IndexCartesian, A::ApplyArray{…}, kr::InfiniteArrays.OneToInf{…}, jr::Int64)
    @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:414
  [7] _getindex
    @ .\multidimensional.jl:889 [inlined]
  [8] getindex
    @ .\abstractarray.jl:1291 [inlined]
  [9] map (repeats 2 times)
    @ .\tuple.jl:342 [inlined]
 [10] _mul_getindex(args::Tuple{ApplyArray{…}, ApplyArray{…}}, k::Int64, j::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:218
 [11] getindex(::ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{…}, ApplyArray{…}}}, ::Int64, ::Int64)
    @ LazyArrays C:\Users\User\.julia\packages\LazyArrays\MLFsy\src\linalg\mul.jl:184
 [12] top-level scope
    @ REPL[47]:1
Some type information was truncated. Use `show(err)` to see complete types.

I must be at the limit of my understanding of the layout system since I'm at a bit of a loss here. I'm trying to get this fixed since it's probably related to all the issues in #109.

Also not sure why A meets the ambiguity but AA doesn't (e: ah, because AA hits some bidiagonal layout methods instead - A is not getting the correct rowsupport instead for example

julia> LazyArrays._mul_args_rowsupport(A.args[1], k)
OneToInf()

julia> LazyArrays._mul_args_rowsupport(AA.args[1], k)
1:2
DanielVandH commented 1 month ago

The issue with the rowsupport mentioned seems to be an issue with finding the row support of a banded layout:

julia> C = A.args[1].args[2]
ℵ₀×ℵ₀ UpperTriangular{Float64, 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}}}} with indices OneToInf()×OneToInf():
 0.68313  0.323669  0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …
 0.0      0.624188  0.382235  0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0      0.0       0.593394  0.410901  0.0       0.0       0.0       0.0       0.0       0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0      0.0       0.0       0.574777  0.428194  0.0       0.0       0.0       0.0       0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     
 0.0      0.0       0.0       0.0       0.562341  0.439818  0.0       0.0       0.0       0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0      0.0       0.0       0.0       0.0       0.553451  0.448184  0.0       0.0       0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  …
 0.0      0.0       0.0       0.0       0.0       0.0       0.54678   0.454499  0.0       0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0      0.0       0.0       0.0       0.0       0.0       0.0       0.54159   0.459438  0.0       0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0      0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.537436  0.463407  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 ⋮                                                ⋮                                                 ⋮                        ⋮                        ⋱

julia> lay = MemoryLayout(C)
BandedMatrices.BandedLayout

julia> rowsupport(lay, C, 1)
OneToInf()

Probably similar for rowsupport. I'm guessing defining

rowsupport(::BandedLayout, A, k) = ((bl, bu) = bandwidths(A); return (k - bl):(k + bl)) # need extra care for k=1 or k = length(A) i think

could be a good solution, unless I'm missing something obvious here.

-

Oh but it already has these methods..

https://github.com/JuliaLinearAlgebra/BandedMatrices.jl/blob/3401f05e65eacd1f97cc8fa6cd5d8388eb9bd3e5/src/generic/AbstractBandedMatrix.jl#L236-L237

just have to make it get to those methods I guess

DanielVandH commented 1 month ago

... It's because I did

ArrayLayouts.MemoryLayout(::Type{<:UpperTriangular{V, <:InfiniteLinearAlgebra.AdaptiveCholeskyFactors{T, <:BandedMatrices.AbstractBandedMatrix}}}) where {V, T} = BandedMatrices.BandedLayout

instead of

ArrayLayouts.MemoryLayout(::Type{<:UpperTriangular{V, <:InfiniteLinearAlgebra.AdaptiveCholeskyFactors{T, <:BandedMatrices.AbstractBandedMatrix}}}) where {V, T} = BandedMatrices.BandedLayout()

(notice the () at the end). Yikes. Works with this:

julia> A
((ℵ₀×ℵ₀ Diagonal{Float64, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, 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}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) * ((ℵ₀×ℵ₀ Diagonal{Float64, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{Bool, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
 1.0  -0.219853  -0.343287    ⋅          ⋅          ⋅           ⋅           ⋅           ⋅           ⋅           ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    0.829116  -0.170677  -0.421318    ⋅          ⋅           ⋅           ⋅           ⋅           ⋅           ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅     ⋅         0.762775  -0.134667  -0.461939    ⋅           ⋅           ⋅           ⋅           ⋅           ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅     ⋅          ⋅         0.726547  -0.110501  -0.487195     ⋅           ⋅           ⋅           ⋅           ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅     ⋅          ⋅          ⋅         0.703486  -0.0935081  -0.50448      ⋅           ⋅           ⋅           ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅     ⋅          ⋅          ⋅          ⋅         0.687467   -0.0809808  -0.51707      ⋅           ⋅           ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …
  ⋅     ⋅          ⋅          ⋅          ⋅          ⋅          0.675677   -0.0713855  -0.526654     ⋅           ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅     ⋅          ⋅          ⋅          ⋅          ⋅           ⋅          0.666631   -0.0638094  -0.534195     ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅     ⋅          ⋅          ⋅          ⋅          ⋅           ⋅           ⋅          0.659469   -0.0576796  -0.540286   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
 ⋮                                                 ⋮                                                           ⋮                             ⋮         ⋱

Will try and get some PRs in for this, although its my MemoryLayout that makes those rowsupports go incorrectly in the first place so more investigating is to be done, since the original issue still wouldn't work:

julia> inv(A)
inv(((ℵ₀×ℵ₀ Diagonal{Float64, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, 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}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) * ((ℵ₀×ℵ₀ Diagonal{Float64, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{Bool, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
Error showing value of type ApplyArray{Float64, 2, typeof(inv), Tuple{ApplyArray{Float64, 2, typeof(*), Tuple{ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}, UpperTriangular{Float64, 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}}}}}}, ApplyArray{Float64, 2, typeof(*), Tuple{Diagonal{Float64, Fill{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}, UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, BroadcastVector{Float64, typeof(-), Tuple{Bool, SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}, BroadcastVector{Float64, typeof(-), Tuple{SubArray{Float64, 1, ClassicalOrthogonalPolynomials.CholeskyJacobiData{Float64}, Tuple{Base.Slice{InfiniteArrays.OneToInf{Int64}}, Int64}, false}}}}}}}}}}}}:
ERROR: StackOverflowError:
SYSTEM (REPL): showing an error caused an error
ERROR: InterruptException:
Stacktrace:
  [1] unsafe_write(to::IOBuffer, p::Ptr{UInt8}, nb::UInt64)
    @ Base .\iobuffer.jl:457
  [2] unsafe_write
    @ .\io.jl:431 [inlined]
  [3] write
    @ .\strings\io.jl:248 [inlined]
  [4] print
    @ .\strings\io.jl:250 [inlined]
  [5] show_typeparams(io::IOContext{IOBuffer}, env::Core.SimpleVector, orig::Core.SimpleVector, wheres::Vector{TypeVar})
    @ Base .\show.jl:696
  [6] show_datatype(io::IOContext{IOBuffer}, x::DataType, wheres::Vector{TypeVar})
    @ Base .\show.jl:1130
  [7] show_datatype
    @ .\show.jl:1077 [inlined]
  [8] _show_type(io::IOContext{IOBuffer}, x::Type)
    @ Base .\show.jl:961
  [9] show(io::IOContext{IOBuffer}, x::Type)
    @ Base .\show.jl:953
 [10] show_typeparams(io::IOContext{IOBuffer}, env::Core.SimpleVector, orig::Core.SimpleVector, wheres::Vector{TypeVar})
    @ Base .\show.jl:710
 [11] show_datatype(io::IOContext{IOBuffer}, x::DataType, wheres::Vector{TypeVar})
    @ Base .\show.jl:1130
 [12] show_datatype
    @ .\show.jl:1077 [inlined]
 [13] _show_type(io::IOContext{IOBuffer}, x::Type)
    @ Base .\show.jl:961
 [14] show(io::IOContext{IOBuffer}, x::Type)
    @ Base .\show.jl:953
 [15] sprint(f::Function, args::Type; context::IOContext{IOBuffer}, sizehint::Int64)
    @ Base .\strings\io.jl:112
 [16] sprint
    @ .\strings\io.jl:107 [inlined]
 [17] #print_type_bicolor#564
    @ .\show.jl:2672 [inlined]
 [18] show_tuple_as_call(out::IOBuffer, name::Symbol, sig::Type; demangle::Bool, kwargs::Base.Iterators.Zip{Tuple{Vector{Symbol}, Vector{Any}}}, argnames::Vector{Symbol}, qualified::Bool, hasfirst::Bool)
    @ Base .\show.jl:2539
 [19] show_tuple_as_call
    @ .\show.jl:2507 [inlined]
 [20] show_spec_sig(io::IOBuffer, m::Method, sig::Type)
    @ Base.StackTraces .\stacktraces.jl:260
 [21] show_spec_linfo(io::IOBuffer, frame::Base.StackTraces.StackFrame)
    @ Base.StackTraces .\stacktraces.jl:232
 [22] show(io::IOBuffer, frame::Base.StackTraces.StackFrame)
    @ Base.StackTraces .\stacktraces.jl:270
 [23] sprint(f::Function, args::Base.StackTraces.StackFrame; context::Nothing, sizehint::Int64)
    @ Base .\strings\io.jl:114
 [24] sprint
    @ .\strings\io.jl:107 [inlined]
 [25] _collapse_repeated_frames(trace::Vector{Any})
    @ Base .\errorshow.jl:863
 [26] process_backtrace(t::Vector{Base.StackTraces.StackFrame}, limit::Int64; skipC::Bool)
    @ Base .\errorshow.jl:961
 [27] process_backtrace (repeats 2 times)
    @ .\errorshow.jl:908 [inlined]
 [28] show_backtrace(io::IOContext{Base.TTY}, t::Vector{Base.StackTraces.StackFrame})
    @ Base .\errorshow.jl:784
 [29] showerror(io::IOContext{Base.TTY}, ex::StackOverflowError, bt::Vector{Base.StackTraces.StackFrame}; backtrace::Bool)
    @ Base .\errorshow.jl:97
 [30] show_exception_stack(io::IOContext{Base.TTY}, stack::Base.ExceptionStack)
    @ Base .\errorshow.jl:975
 [31] display_error(io::IOContext{Base.TTY}, stack::Base.ExceptionStack)
    @ Base .\client.jl:111
 [32] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [33] invokelatest
    @ .\essentials.jl:889 [inlined]
 [34] repl_display_error(errio::IO, errval::Any)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:293
 [35] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:310
 [36] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:284
 [37] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [38] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:282
 [39] (::REPL.var"#do_respond#80"{Bool, Bool, REPL.var"#93#103"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:911
 [40] (::VSCodeServer.var"#103#106"{REPL.var"#do_respond#80"{Bool, Bool, REPL.var"#93#103"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt}})(mi::REPL.LineEdit.MIState, buf::IOBuffer, ok::Bool)
    @ VSCodeServer c:\Users\User\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:122
 [41] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [42] invokelatest
    @ .\essentials.jl:889 [inlined]
 [43] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\LineEdit.jl:2656
 [44] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:1312
 [45] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:386
DanielVandH commented 1 month ago

Ignoring all the above where I keep tripping over myself, the issue is reducible to something like

using ClassicalOrthogonalPolynomials, SemiclassicalOrthogonalPolynomials, ArrayLayouts
A = SemiclassicalJacobi(2.0, 2.0, 1.0, 1.0) \ SemiclassicalJacobi(2.0, 1.0, 0.0, 1.0)
B = A.args[1].args[2]
julia> B = A.args[1].args[2]
ℵ₀×ℵ₀ LinearAlgebra.UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrices.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}}}} with indices OneToInf()×OneToInf():
 0.68313  0.323669   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …
  ⋅       0.624188  0.382235   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅        0.593394  0.410901   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅        ⋅         ⋅        0.574777  0.428194   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅         ⋅         ⋅        0.562341  0.439818   ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅         ⋅         ⋅         ⋅        0.553451  0.448184   ⋅         ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …
  ⋅        ⋅         ⋅         ⋅         ⋅         ⋅        0.54678   0.454499   ⋅         ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        0.54159   0.459438   ⋅         ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅        ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        0.537436  0.463407   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
 ⋮                                                ⋮                                                 ⋮                        ⋮                        ⋱
julia> inv(B)
ERROR: MethodError: no method matching Matrix{…}(::LinearAlgebra.UniformScaling{…}, ::Tuple{…})

Closest candidates are:
  Array{T, N}(::Missing, ::Any...) where {T, N}
   @ Base baseext.jl:43
  Array{T, N}(::Nothing, ::Any...) where {T, N}
   @ Base baseext.jl:42
  Matrix{T}(::UndefInitializer, ::Tuple{Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}, Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}}) where T
   @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:7
  ...

Stacktrace:
 [1] inv(A::LinearAlgebra.UpperTriangular{Float64, InfiniteLinearAlgebra.AdaptiveCholeskyFactors{…}})
   @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:797
 [2] top-level scope
   @ c:\Users\User\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\test\runtests.jl:619
Some type information was truncated. Use `show(err)` to see complete types.

since inv(B) is trying to use this method

https://github.com/JuliaLang/julia/blob/319d8900dbb2566eef9f7e594e7a1e875dfe819a/stdlib/LinearAlgebra/src/triangular.jl#L1020-L1028

for whatever reason. Probably related to this example

julia> using InfiniteArrays, BandedMatrices, LinearAlgebra
julia> B = BandedMatrices._BandedMatrix(Ones(∞)', ∞, 0, 0) |> UpperTriangular
ℵ₀×ℵ₀ UpperTriangular{Float64, BandedMatrix{Float64, Adjoint{Float64, Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}, InfiniteArrays.OneToInf{Int64}}} with indices OneToInf()×OneToInf():
 1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
 ⋮                        ⋮                        ⋮                        ⋮                        ⋮                        ⋮                        ⋱

julia> inv(B)
ERROR: MethodError: no method matching Matrix{Float64}(::UniformScaling{Bool}, ::Tuple{Infinities.InfiniteCardinal{0}, Infinities.InfiniteCardinal{0}})

Closest candidates are:
  Array{T, N}(::Missing, ::Any...) where {T, N}
   @ Base baseext.jl:43
  Array{T, N}(::Nothing, ::Any...) where {T, N}
   @ Base baseext.jl:42
  Matrix{T}(::UndefInitializer, ::Tuple{Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}, Union{Infinities.InfiniteCardinal{0}, Infinities.Infinity}}) where T
   @ InfiniteArrays C:\Users\User\.julia\packages\InfiniteArrays\7v9AI\src\infarrays.jl:7
  ...

Stacktrace:
 [1] inv(A::UpperTriangular{Float64, BandedMatrix{Float64, Adjoint{Float64, Ones{…}}, InfiniteArrays.OneToInf{Int64}}})
   @ LinearAlgebra C:\Users\User\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\triangular.jl:797
 [2] top-level scope
   @ REPL[20]:1
Some type information was truncated. Use `show(err)` to see complete types.
TSGut commented 1 month ago

We can always overload inv to hit a method we custom make for the connection coefficients since they have their own structs underneath the UpperTriangular and are thus easily targeted. That'll make it much easier to end up with exactly the layout you want.

That method would belong in ClassicalOPs though.

DanielVandH commented 6 days ago

This works now