search

JuliaApproximation / SemiclassicalOrthogonalPolynomials.jl

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

Slicing with infinite range tries to evaluate entire matrix #120

Open DanielVandH opened 1 month ago

DanielVandH commented 1 month ago

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.

TSGut commented 1 month ago

Looking at the output for the case that works it's probably just overloading views anyway somewhere in InfiniteArrays or InfiniteLinearAlgebra. You could use Debugger to find where that happens and add a line to that package for appropriate form ApplyLayouts as well.