JuliaSIMD / LoopVectorization.jl

Macro(s) for vectorizing loops.
MIT License
742 stars 66 forks source link

Incorrect results using @turbo with linear array indexing #512

Open ajshephard opened 12 months ago

ajshephard commented 12 months ago

There was some discussion about this issue here. Essentially, when using linear indexing of an array @turbo appears to start indexing the array incorrectly when a offset term is used. An example of this (code and commentary copied directly from the discourse discussion) is below.

julia> H = collect(reshape(1.0:6.0, 3, 2))
3×2 Matrix{Float64}:
 1.0  4.0
 2.0  5.0
 3.0  6.0

julia> su = zeros(3);

julia> for k = 1:1
           ind = 0 # LinearIndices((3, 2))[1, 1] - 1
           for ij in 1:3
               su[ij] += H[ij + ind]
           end
       end

julia> su
3-element Vector{Float64}:
 1.0
 2.0
 3.0

julia> su = zeros(3);

julia> for k = 1:1
           ind = 0
           @turbo for ij in 1:3
               su[ij] += H[ij + ind]
           end
       end

julia> su
3-element Vector{Float64}:
 6.94274366677856e-310
 1.0
 2.0

so H is being indexed at 0, 1, 2 instead of 1, 2, 3. The effect doesn’t seem to happen on the left hand expression e.g. su[ij + ind] and doesn’t seem to shift more with more additions e.g. H[ij + ind + ind]. It does it correctly if the index was added to a literal 0 instead of ind or if not added:

julia> su = zeros(3);

julia> for k = 1:1
           @turbo for ij in 1:3
               su[ij] += H[ij + 0]
           end
       end

julia> su
3-element Vector{Float64}:
 1.0
 2.0
 3.0

julia> su = zeros(3);

julia> for k = 1:1
           @turbo for ij in 1:3
               su[ij] += H[ij]
           end
       end

julia> su
3-element Vector{Float64}:
 1.0
 2.0
 3.0
chriselrod commented 12 months ago

Probably a bug in some of the offsetting code. If anyone wants to take a look at it, I'd be happy to answer any questions.

I am unlikely to get around to fixing this myself.