JuliaSIMD / StrideArrays.jl

Library supporting the ArrayInterface.jl strided array interface.
MIT License
54 stars 9 forks source link

Error broadcasting a stride array with a unitrange: #76

Open MasonProtter opened 1 year ago

MasonProtter commented 1 year ago

This works fine with StrideArraysCore.jl, but not with StrideArrays.jl

#+begin_src julia
using StrideArrays

let A = StrideArray{Float64}(undef, 10, 10)
    v = 1:10
    A .= v .* v'
end
#+end_src

#+RESULTS:
:RESULTS:
# [goto error]
#+begin_example
conversion to pointer not defined for UnitRange{Int64}

Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] unsafe_convert(#unused#::Type{Ptr{Int64}}, a::UnitRange{Int64})
    @ Base ./pointer.jl:67
  [3] unsafe_convert(#unused#::Type{Ptr{Int64}}, A::Adjoint{Int64, UnitRange{Int64}})
    @ LinearAlgebra ~/julia/usr/share/julia/stdlib/v1.9/LinearAlgebra/src/adjtrans.jl:321
  [4] pointer(x::Adjoint{Int64, UnitRange{Int64}})
    @ Base ./abstractarray.jl:1240
  [5] unsafe_convert(#unused#::Type{Ptr{Int64}}, A::LoopVectorization.LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{Int64}}})
    @ LoopVectorization ~/.julia/packages/LoopVectorization/5ukqQ/src/broadcast.jl:63
  [6] pointer(x::LoopVectorization.LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{Int64}}})
    @ Base ./abstractarray.jl:1240
  [7] memory_reference
    @ ~/.julia/packages/LayoutPointers/qGkBo/src/stridedpointers.jl:21 [inlined]
  [8] memory_reference
    @ ~/.julia/packages/LayoutPointers/qGkBo/src/stridedpointers.jl:18 [inlined]
  [9] memory_reference
    @ ~/.julia/packages/LoopVectorization/5ukqQ/src/broadcast.jl:197 [inlined]
 [10] stridedpointer_preserve
    @ ~/.julia/packages/LayoutPointers/qGkBo/src/stridedpointers.jl:100 [inlined]
 [11] macro expansion
    @ ~/.julia/packages/LoopVectorization/5ukqQ/src/broadcast.jl:671 [inlined]
 [12] vmaterialize!
    @ ~/.julia/packages/LoopVectorization/5ukqQ/src/broadcast.jl:671 [inlined]
 [13] vmaterialize!
    @ ~/.julia/packages/LoopVectorization/5ukqQ/src/broadcast.jl:757 [inlined]
 [14] _materialize!
    @ ~/.julia/packages/StrideArrays/lTVQX/src/broadcast.jl:205 [inlined]
 [15] materialize!(dest::StrideArray{Float64, 2, (1, 2), Tuple{Int64, Int64}, Tuple{Nothing, Nothing}, Tuple{StaticInt{1}, StaticInt{1}}, Vector{Float64}}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(*), Tuple{UnitRange{Int64}, Adjoint{Int64, UnitRange{Int64}}}})
    @ StrideArrays ~/.julia/packages/StrideArrays/lTVQX/src/broadcast.jl:212
 [16] top-level scope
    @ In[37]:5
#+end_example
:END:
chriselrod commented 1 year ago

StrideArrays.jl pirates broadcasting + StrideArrays, using LoopVectorization.

The specific problem here is the adjoint unit range + loop vectorization's broadcast handling.

julia> using LoopVectorization

julia> B = rand(10,10); v = 1:10;

julia> @turbo B .= v .* v'
ERROR: conversion to pointer not defined for UnitRange{Int64}
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] unsafe_convert(::Type{Ptr{Int64}}, a::UnitRange{Int64})
    @ Base ./pointer.jl:67
  [3] unsafe_convert(::Type{Ptr{Int64}}, A::Adjoint{Int64, UnitRange{Int64}})
    @ LinearAlgebra ~/Documents/languages/juliamaster/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/adjtrans.jl:356
  [4] pointer(x::Adjoint{Int64, UnitRange{Int64}})
    @ Base ./abstractarray.jl:1232
  [5] unsafe_convert(::Type{Ptr{Int64}}, A::LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{Int64}}})
    @ LoopVectorization ~/.julia/dev/LoopVectorization/src/broadcast.jl:63
  [6] pointer(x::LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{Int64}}})
    @ Base ./abstractarray.jl:1232
  [7] memory_reference(::ArrayInterface.CPUPointer, A::LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{Int64}}})
    @ LayoutPointers ~/.julia/dev/LayoutPointers/src/stridedpointers.jl:21 [inlined]
  [8] memory_reference(A::LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{Int64}}})
    @ LayoutPointers ~/.julia/dev/LayoutPointers/src/stridedpointers.jl:18 [inlined]
  [9] memory_reference(fb::LoopVectorization.ForBroadcast{Int64, 2, LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{Int64}}}})
    @ LoopVectorization ~/.julia/dev/LoopVectorization/src/broadcast.jl:197 [inlined]
 [10] stridedpointer_preserve(A::LoopVectorization.ForBroadcast{Int64, 2, LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{…}}}})

while without the adjoint, it works

julia> b = rand(10);

julia> @turbo b .= v .* v
10-element Vector{Float64}:
   1.0
   4.0
   9.0
  16.0
  25.0
  36.0
  49.0
  64.0
  81.0
 100.0

Note from the stack trace it has been converted into a "low-dim-array":

 [10] stridedpointer_preserve(A::LoopVectorization.ForBroadcast{Int64, 2, LowDimArray{(false, true), Int64, 2, Adjoint{Int64, UnitRange{…}}}})
    @ LayoutPointers ~/.julia/dev/LayoutPointers/src/stridedpointers.jl:100 [inlined]

Although we also get that for UnitRanges:

julia> LoopVectorization.stridedpointer_preserve(v')
ERROR: "Memory access for Adjoint{Int64, UnitRange{Int64}} not implemented yet."
Stacktrace:
 [1] memory_reference(::ArrayInterface.CPUIndex, A::Adjoint{Int64, UnitRange{Int64}})
   @ LayoutPointers ~/.julia/dev/LayoutPointers/src/stridedpointers.jl:62 [inlined]
 [2] memory_reference(A::Adjoint{Int64, UnitRange{Int64}})
   @ LayoutPointers ~/.julia/dev/LayoutPointers/src/stridedpointers.jl:18 [inlined]
 [3] stridedpointer_preserve(A::Adjoint{Int64, UnitRange{Int64}})
   @ LayoutPointers ~/.julia/dev/LayoutPointers/src/stridedpointers.jl:100
 [4] top-level scope
   @ REPL[18]:1

Probably a few routes you could take to fixing this. Two ideas:

  1. Add a transposed FastRange type (i.e. what stridedpointer(1:10) returns, and overloads to produce it here.
  2. Explicit support within the add_broadcast function: https://github.com/JuliaSIMD/LoopVectorization.jl/blob/5a4d58f66aef373371d6e0f45ecef60589952ac7/src/broadcast.jl#L480-L508