Slicing with infinite range tries to evaluate entire matrix #120
Description
I don't know what repo this issue should actually belong to, but see:
julia> using SemiclassicalOrthogonalPolynomials
julia> P = SemiclassicalJacobi(2.0, -1/2, -1.0, -1/2);
julia> P0 = SemiclassicalJacobi(2.0, -1/2, 0.0, -1/2);
julia> R = P0 \ P
(ℵ₀×ℵ₀ LinearAlgebra.UpperTriangular{Float64, LinearAlgebra.Diagonal{Float64, FillArrays.Ones{Float64, 1, Tuple{InfiniteArrays.OneToInf{Int64}}}}} with indices OneToInf()×OneToInf()) * (ℵ₀×ℵ₀ LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.ApplyArray{Float64, 1, typeof(vcat), Tuple{Float64, SubArray{Float64, 1, LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}}}, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, OrthogonalPolynomialRatio{Float64, SemiclassicalJacobi{Float64}}}}}}}}}, Tuple{InfiniteLinearAlgebra.InfBandCartesianIndices}, false}}}, SubArray{Float64, 1, LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}}}, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{LazyArrays.BroadcastVector{Float64, typeof(*), Tuple{Float64, LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, OrthogonalPolynomialRatio{Float64, SemiclassicalJacobi{Float64}}}}}}}}}, Tuple{InfiniteLinearAlgebra.InfBandCartesianIndices}, false}} with indices OneToInf()×OneToInf()) with indices OneToInf()×OneToInf():
1.0 0.63662 ⋅ ⋅ ⋅ ⋅ …
⋅ -0.307758 0.461696 ⋅ ⋅ ⋅
⋅ ⋅ -0.351058 0.437497 ⋅ ⋅
⋅ ⋅ ⋅ -0.36663 0.426727 ⋅
⋅ ⋅ ⋅ ⋅ -0.374557 0.420662
⋅ ⋅ ⋅ ⋅ ⋅ -0.379357 …
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋮ ⋮ ⋱
julia> R[2:3, 2:3] # ok
2×2 Matrix{Float64}:
-0.307758 0.461696
0.0 -0.351058
julia> R[2:end, 2:end]
ERROR: InterruptException:
Stacktrace:
[1] checkindex(::Type{Bool}, inds::InfiniteArrays.OneToInf{Int64}, i::Int64)
@ Base .\abstractarray.jl:761
[2] checkbounds
@ .\abstractarray.jl:687 [inlined]
[3] _vec_resizedata!(B::LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}, n::Int64)
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:232
[4] resizedata!
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:251 [inlined]
[5] resizedata!
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:219 [inlined]
[6] getindex(A::LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{…}, AbstractVector{…}}, I::CartesianIndex{1})
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:206
[7] _broadcast_getindex
@ .\broadcast.jl:662 [inlined]
[8] _getindex
@ .\broadcast.jl:705 [inlined]
[9] _broadcast_getindex
@ .\broadcast.jl:681 [inlined]
[10] _getindex
@ .\broadcast.jl:705 [inlined]
[11] _broadcast_getindex(bc::Base.Broadcast.Broadcasted{…}, I::Int64)
@ Base.Broadcast .\broadcast.jl:681
[12] getindex
@ .\broadcast.jl:636 [inlined]
[13] getindex
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\lazybroadcasting.jl:93 [inlined]
[14] getindex
@ .\subarray.jl:290 [inlined]
[15] _broadcast_getindex
@ .\broadcast.jl:675 [inlined]
[16] _getindex(args::Tuple{Base.Broadcast.Extruded{…}, Base.Broadcast.Extruded{…}}, I::Int64)
@ Base.Broadcast .\broadcast.jl:705
[17] _broadcast_getindex(bc::Base.Broadcast.Broadcasted{Nothing, Tuple{…}, typeof(*), Tuple{…}}, I::Int64)
@ Base.Broadcast .\broadcast.jl:681
[18] getindex
@ .\broadcast.jl:636 [inlined]
[19] macro expansion
@ .\broadcast.jl:1004 [inlined]
[20] macro expansion
@ .\simdloop.jl:77 [inlined]
[21] copyto!
@ .\broadcast.jl:1003 [inlined]
[22] copyto!
@ .\broadcast.jl:956 [inlined]
[23] copy
@ .\broadcast.jl:928 [inlined]
[24] materialize
@ .\broadcast.jl:903 [inlined]
[25] broadcast(::typeof(*), ::SubArray{Float64, 1, LazyArrays.BroadcastVector{…}, Tuple{…}, false}, ::Vector{Float64})
@ Base.Broadcast .\broadcast.jl:841
[26] layout_broadcasted(::LazyArrays.BroadcastLayout{…}, ::LazyArrays.PaddedColumns{…}, ::typeof(*), A::LazyArrays.BroadcastVector{…}, B::LazyArrays.ApplyArray{…})
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\padded.jl:328
[27] layout_broadcasted(op::Function, A::LazyArrays.BroadcastVector{…}, B::LazyArrays.ApplyArray{…})
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:370
[28] broadcasted
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\lazyconcat.jl:556 [inlined]
[29] broadcasted
@ .\broadcast.jl:1347 [inlined]
[30] copy(M::ArrayLayouts.Lmul{…})
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\diagonal.jl:19
[31] copy
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\padded.jl:508 [inlined]
[32] materialize
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:137 [inlined]
[33] mul
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:138 [inlined]
[34] *
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:179 [inlined]
[35] _tail_simplify
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\linalg\mul.jl:352 [inlined]
[36] _most_simplify
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\linalg\mul.jl:351 [inlined]
[37] _simplify(::typeof(*), ::LinearAlgebra.Diagonal{…}, ::ClassicalOrthogonalPolynomials.Clenshaw{…}, ::SubArray{…})
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\linalg\mul.jl:349
[38] simplify
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\linalg\mul.jl:345 [inlined]
[39] simplify
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\linalg\mul.jl:355 [inlined]
[40] copy
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\linalg\mul.jl:360 [inlined]
[41] copy
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\ext\LazyArraysBandedMatricesExt.jl:638 [inlined]
[42] materialize
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:137 [inlined]
[43] mul
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:138 [inlined]
[44] dot(x::SubArray{…}, A::LazyArrays.ApplyArray{…}, y::SubArray{…})
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:474
[45] dot
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:475 [inlined]
[46] lanczos!(Ns::UnitRange{…}, X::LazyBandedMatrices.SymTridiagonal{…}, W::LazyArrays.ApplyArray{…}, γ::LazyArrays.CachedArray{…}, β::LazyArrays.CachedArray{…}, R::LinearAlgebra.UpperTriangular{…})
@ ClassicalOrthogonalPolynomials C:\Users\danjv\.julia\packages\ClassicalOrthogonalPolynomials\3YYrU\src\lanczos.jl:13
[47] resizedata!(L::ClassicalOrthogonalPolynomials.LanczosData{Float64}, N::Int64)
@ ClassicalOrthogonalPolynomials C:\Users\danjv\.julia\packages\ClassicalOrthogonalPolynomials\3YYrU\src\lanczos.jl:62
[48] resizedata!
@ C:\Users\danjv\.julia\packages\ClassicalOrthogonalPolynomials\3YYrU\src\lanczos.jl:132 [inlined]
[49] _lanczos_getindex
@ C:\Users\danjv\.julia\packages\ClassicalOrthogonalPolynomials\3YYrU\src\lanczos.jl:137 [inlined]
[50] getindex
@ C:\Users\danjv\.julia\packages\ClassicalOrthogonalPolynomials\3YYrU\src\lanczos.jl:145 [inlined]
[51] getindex
@ .\subarray.jl:290 [inlined]
[52] _broadcast_getindex
@ .\broadcast.jl:675 [inlined]
[53] _getindex (repeats 4 times)
@ .\broadcast.jl:705 [inlined]
[54] _broadcast_getindex
@ .\broadcast.jl:681 [inlined]
--- the last 2 lines are repeated 1 more time ---
[57] _getindex
@ .\broadcast.jl:706 [inlined]
[58] _broadcast_getindex
@ .\broadcast.jl:681 [inlined]
[59] getindex
@ .\broadcast.jl:636 [inlined]
[60] macro expansion
@ .\broadcast.jl:1004 [inlined]
[61] macro expansion
@ .\simdloop.jl:77 [inlined]
[62] copyto!
@ .\broadcast.jl:1003 [inlined]
[63] copyto!
@ .\broadcast.jl:956 [inlined]
[64] materialize!
@ .\broadcast.jl:914 [inlined]
[65] materialize!(dest::SubArray{…}, bc::Base.Broadcast.Broadcasted{…})
@ Base.Broadcast .\broadcast.jl:911
[66] copyto!_layout(::ArrayLayouts.DenseColumnMajor, ::LazyArrays.BroadcastLayout{…}, dest::SubArray{…}, bc::SubArray{…})
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\lazybroadcasting.jl:24
[67] copyto!_layout
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\ArrayLayouts.jl:267 [inlined]
[68] copyto!(dest::SubArray{…}, src::SubArray{…})
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\ArrayLayouts.jl:276
[69] vcat_copyto!(::SubArray{Float64, 1, Vector{…}, Tuple{…}, true}, ::FillArrays.Fill{Float64, 1, Tuple{…}}, ::Vararg{Any})
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\lazyconcat.jl:257
[70] copyto!_layout
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\lazyconcat.jl:226 [inlined]
[71] copyto!_layout
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\ArrayLayouts.jl:267 [inlined]
[72] copyto!(dest::SubArray{…}, src::SubArray{…})
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\ArrayLayouts.jl:276
[73] cache_filldata!(K::LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{…}, AbstractVector{…}}, inds::UnitRange{Int64})
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\lazyoperations.jl:442
[74] _vec_resizedata!(B::LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}, n::Int64)
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:243
[75] resizedata!
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:251 [inlined]
[76] resizedata!
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:219 [inlined]
[77] getindex(A::LazyArrays.Accumulate{Float64, 1, typeof(*), Vector{…}, AbstractVector{…}}, I::CartesianIndex{1})
@ LazyArrays C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\cache.jl:206
[78] _broadcast_getindex
@ .\broadcast.jl:662 [inlined]
[79] _getindex
@ .\broadcast.jl:706 [inlined]
[80] _getindex
@ .\broadcast.jl:705 [inlined]
[81] _broadcast_getindex
@ .\broadcast.jl:681 [inlined]
[82] _getindex
@ .\broadcast.jl:706 [inlined]
[83] _broadcast_getindex
@ .\broadcast.jl:681 [inlined]
[84] getindex(bc::Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{…}, Tuple{…}, typeof(-), Tuple{…}}, I::Int64)
@ Base.Broadcast .\broadcast.jl:636
[85] getindex
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\lazybroadcasting.jl:93 [inlined]
[86] getindex(A::LazyBandedMatrices.Bidiagonal{…}, i::Int64, j::Int64)
@ LazyBandedMatrices C:\Users\danjv\.julia\packages\LazyBandedMatrices\dXFwo\src\bidiag.jl:116
[87] _getindex
@ .\abstractarray.jl:1341 [inlined]
[88] getindex
@ .\abstractarray.jl:1291 [inlined]
[89] getindex
@ .\subarray.jl:290 [inlined]
[90] getindex(A::LazyBandedMatrices.Bidiagonal{Float64, LazyArrays.ApplyArray{…}, SubArray{…}}, i::Int64, j::Int64)
@ LazyBandedMatrices C:\Users\danjv\.julia\packages\LazyBandedMatrices\dXFwo\src\bidiag.jl:118
[91] getindex
@ .\subarray.jl:290 [inlined]
[92] macro expansion
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\muladd.jl:167 [inlined]
[93] macro expansion
@ .\simdloop.jl:77 [inlined]
[94] _default_blasmul_loop!
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\muladd.jl:166 [inlined]
[95] default_blasmul!(α::Bool, A::SubArray{…}, B::SubArray{…}, β::Bool, C::LazyArrays.CachedArray{…})
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\muladd.jl:192
[96] materialize!(M::ArrayLayouts.MulAdd{…})
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\muladd.jl:264
[97] muladd!(α::Bool, A::SubArray{…}, B::SubArray{…}, β::Bool, C::LazyArrays.CachedArray{…}; kw::@Kwargs{})
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\muladd.jl:75
[98] mul!(dest::LazyArrays.CachedArray{…}, A::SubArray{…}, B::SubArray{…}, α::Bool, β::Bool)
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:144
[99] mul!(C::LazyArrays.CachedArray{…}, A::SubArray{…}, B::SubArray{…}, α::Bool, β::Bool)
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\mul.jl:393
[100] mul!
@ C:\Users\danjv\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:237 [inlined]
[101] *(A::SubArray{…}, B::SubArray{…})
@ LinearAlgebra C:\Users\danjv\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\matmul.jl:106
[102] sub_materialize
@ C:\Users\danjv\.julia\packages\LazyArrays\6UPSU\src\linalg\mul.jl:261 [inlined]
[103] sub_materialize
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\ArrayLayouts.jl:132 [inlined]
[104] layout_getindex
@ C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\ArrayLayouts.jl:138 [inlined]
[105] getindex(A::LazyArrays.ApplyArray{…}, kr::InfiniteArrays.InfUnitRange{…}, jr::InfiniteArrays.InfUnitRange{…})
@ ArrayLayouts C:\Users\danjv\.julia\packages\ArrayLayouts\31idh\src\ArrayLayouts.jl:153
Some type information was truncated. Use `show(err)` to see complete types.The [2:end, 2:end] slice is trying to evaluating all the entries. Should this work properly? It's handled fine for for e.g.
julia> P = SemiclassicalJacobi(2.0, -1/2, 1.0, -1/2);
julia> P0 = SemiclassicalJacobi(2.0, -1/2, 2.0, -1/2);
julia> R = P0 \ P;
julia> R[2:end, 2:end]
ℵ₀×ℵ₀ view(::LinearAlgebra.UpperTriangular{Float64, LazyArrays.BroadcastMatrix{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.AdaptiveCholeskyFactors{Float64, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, LazyBandedMatrices.SymTridiagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{Bool, ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}, LazyArrays.BroadcastVector{Float64, typeof(-), Tuple{ClassicalOrthogonalPolynomials.LanczosJacobiBand{Float64}}}}}, Float64}}}, 2:∞, 2:∞) with eltype Float64 with indices OneToInf()×OneToInf():
0.769448 -0.399213 0.0 0.0 0.0 …
0.0 0.701621 -0.445969 0.0 0.0
0.0 0.0 0.667749 -0.471919 0.0
0.0 0.0 0.0 0.647324 -0.488489
0.0 0.0 0.0 0.0 0.633644
0.0 0.0 0.0 0.0 0.0 …
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 …
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
⋮ ⋱The difference in this case is that, for the latter example,
julia> MemoryLayout(R)
TriangularLayout{'U', 'N', LazyArrays.LazyLayout}()but for the first example
julia> MemoryLayout(R)
LazyArrays.ApplyLayout{typeof(*)}()(And their types are different.) Any idea what the issue is here? I can get around this using @view for now but seems like a weird issue that it's not automatically creating the view for this infinite range as I'd expect when not using b=-1 anywhere.