Issue with r'*X*r for X::KronTrav and r::BroadcastVector

[DanielVandH]

Not really 100% sure that this is a KronTrav issue, but I've only ever encountered it with KronTrav so far. Weird issue with r'*X*r whenever X is a KronTrav and r is a BroadcastVector.

julia> using LazyBandedMatrices, InfiniteRandomArrays, LazyArrays, InfiniteLinearAlgebra

julia> A, B = [InfRandSymTridiagonal() for _ in 1:2];

julia> X = KronTrav(A, B);

julia> r = InfiniteLinearAlgebra.pad([1.0], axes(X, 2));

julia> r'*X*r # works fine
0.12917024894291182

julia> r = 1.0 * r;

julia> typeof(r)
BroadcastVector{Float64, typeof(*), Tuple{Float64, ApplyArray{Float64, 1, typeof(setindex), Tuple{Zeros{Float64, 1, Tuple{BlockedOneTo{Int64, RangeCumsum{Int64, OneToInf{Int64}}}}}, Vector{Float64}, OneTo{Int64}}}}} (alias for BroadcastArray{Float64, 1, typeof(*), Tuple{Float64, ApplyArray{Float64, 1, typeof(Base.setindex), Tuple{Zeros{Float64, 1, Tuple{BlockedOneTo{Int64, ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}, Array{Float64, 1}, Base.OneTo{Int64}}}}})

julia> r'*X*r # no longer works fine
ERROR: InterruptException:
Stacktrace:
  [1] MemoryLayout
    @ C:\Users\djv23\.julia\packages\BlockArrays\8VbzB\src\blocklinalg.jl:81 [inlined]
  [2] MemoryLayout
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\memorylayout.jl:169 [inlined]
  [3] colsupport(lay::LazyArrays.PaddedColumns{…}, A::ApplyArray{…}, j::Int64)
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\ltmzk\src\padded.jl:88
  [4] colsupport
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\memorylayout.jl:662 [inlined]
  [5] _broadcast_colsupport
    @ C:\Users\djv23\.julia\packages\LazyArrays\ltmzk\src\lazybroadcasting.jl:199 [inlined]
  [6] _broadcast_getindex_evalf
    @ .\broadcast.jl:678 [inlined]
  [7] _broadcast_getindex
    @ .\broadcast.jl:651 [inlined]
  [8] #17
    @ .\broadcast.jl:1102 [inlined]
  [9] ntuple
    @ .\ntuple.jl:49 [inlined]
 [10] copy
    @ .\broadcast.jl:1102 [inlined]
 [11] materialize
    @ .\broadcast.jl:872 [inlined]
 [12] colsupport(lay::LazyArrays.BroadcastLayout{typeof(*)}, A::BroadcastVector{Float64, typeof(*), Tuple{…}}, j::Int64)
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\ltmzk\src\lazybroadcasting.jl:207
 [13] colsupport
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\memorylayout.jl:662 [inlined]
 [14] _getindex(::Type{…}, M::Mul{…}, ::Tuple{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\mul.jl:40
 [15] getindex
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\mul.jl:76 [inlined]
 [16] _mul_getindex
    @ C:\Users\djv23\.julia\packages\LazyArrays\ltmzk\src\linalg\mul.jl:180 [inlined]
 [17] getindex
    @ C:\Users\djv23\.julia\packages\LazyArrays\ltmzk\src\linalg\mul.jl:184 [inlined]
 [18] dot(x::ApplyArray{Float64, 1, typeof(*), Tuple{…}}, y::BroadcastVector{Float64, typeof(*), Tuple{…}})
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\generic.jl:896
 [19] copy
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\mul.jl:448 [inlined]
 [20] copy
    @ C:\Users\djv23\.julia\packages\LazyArrays\ltmzk\src\linalg\mul.jl:368 [inlined]
 [21] materialize
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\mul.jl:137 [inlined]
 [22] mul
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\mul.jl:138 [inlined]
 [23] *(A::ApplyArray{Float64, 2, typeof(*), Tuple{…}}, B::BroadcastVector{Float64, typeof(*), Tuple{…}})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\B2wRU\src\mul.jl:226
 [24] *(tu::LinearAlgebra.Adjoint{…}, B::KronTrav{…}, v::BroadcastVector{…})
    @ LinearAlgebra C:\Users\djv23\.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\matmul.jl:1117
 [25] top-level scope
    @ REPL[85]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> ApplyArray(*, r', X, r) # this is ok
(((Float64) .* (setindex(ℵ₀-element Zeros{Float64, 1, Tuple{BlockedOneTo{Int64, ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}} with indices BlockedOneTo(ArrayLayouts.RangeCumsum(OneToInf())), 1-element Vector{Float64}, 1-element Base.OneTo{Int64}) with indices BlockedOneTo(ArrayLayouts.RangeCumsum(OneToInf()))) with indices BlockedOneTo(ArrayLayouts.RangeCumsum(OneToInf())))' with indices Base.OneTo(1)×BlockedOneTo(ArrayLayouts.RangeCumsum(OneToInf()))) * (ℵ₀×ℵ₀-blocked ℵ₀×ℵ₀ KronTrav{Float64, 2, Tuple{LinearAlgebra.SymTridiagonal{Float64, InfRandVector{Float64, Type{Float64}, Random.Xoshiro}}, LinearAlgebra.SymTridiagonal{Float64, InfRandVector{Float64, Type{Float64}, Random.Xoshiro}}}, Tuple{BlockedOneTo{Int64, ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}, BlockedOneTo{Int64, ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}}) * ((Float64) .* (setindex(ℵ₀-element Zeros{Float64, 1, Tuple{BlockedOneTo{Int64, ArrayLayouts.RangeCumsum{Int64, InfiniteArrays.OneToInf{Int64}}}}} with indices BlockedOneTo(ArrayLayouts.RangeCumsum(OneToInf())), 1-element Vector{Float64}, 1-element Base.OneTo{Int64}) with indices BlockedOneTo(ArrayLayouts.RangeCumsum(OneToInf()))) with indices BlockedOneTo(ArrayLayouts.RangeCumsum(OneToInf()))):
 0.12917024894291182

julia> r'*X, X*r # these both work
([0.12917024894291182 0.31886141780288035 …  … ], [0.12917024894291182, 0.31886141780288035, 0.09705420877436907, 0.0, 0.23958181444095045, 0.0, 0.0, 0.0, 0.0, 0.0  …  ])

[dlfivefifty]

It looks like the issue is that dot isnot being overloaded. This is a bit surprising since they are both LayoutArray so I would have expected it to be piped through ArrayLayouts.Dot. But that was never fully developed...

[dlfivefifty]

Oh actually two issues: 1.0 * r shouldn't change the type. We would need to look at layout_broadcasted to understand why