JuliaSIMD / JuliaSIMD/LoopVectorization.jl

`@avx` incorrectly handling variable redefinitions between iterations

Offen
#266 5 Kommentare 0 Reaktionen 0 zugewiesene Personen Auf GitHub ansehen
Vorherrschende Sprache
Julia
Sterne
789
Forks
73
PR-Merge-Kennzahlen
Keine gemergten PRs in 30 T.

Beschreibung

I'm not sure if this is something that `@avx` is supposed to be able to handle, but when I write
```julia
using LoopVectorization

function matmul!(u::AbstractVector{T}, A::Tridiagonal{T}, v::AbstractVector{T}) where {T}
@assert length(u) == size(A,1) == size(A,2) == length(v)
dl, d, du = A.dl, A.d, A.du
N = length(u); @assert N > 2

p = zero(T)
c = v[1]
n = v[2]
@inbounds u[1] = d[1] * c + du[1] * n
for i in 2:(N-1)
p = c
c = n
n = v[i+1]
u[i] = dl[i-1] * p + d[i] * c + du[i] * n
end
p = c
c = n
@inbounds u[N] = dl[N-1] * p + d[N] * c
u
end
```
I get correct results compared to `mul!`
```julia
let N = 5
T = Float64
A = Tridiagonal(rand(T, N-1), rand(T, N), rand(T, N-1))
v = rand(T, N)
u1 = Array{T}(undef, N)
u2 = Array{T}(undef, N)

matmul!(u1, A, v)
mul!(u2, A, v)
u1 - u2
end

#+RESULTS:
: 5-element Vector{Float64}:
: 0.0
: 0.0
: 0.0
: 0.0
: 0.0
```

However, sticking `@avx` on the loop makes it get the wrong answer:
```julia
function matmul_avx!(u::AbstractVector{T}, A::Tridiagonal{T}, v::AbstractVector{T}) where {T}
@assert length(u) == size(A,1) == size(A,2) == length(v)
dl, d, du = A.dl, A.d, A.du
N = length(u); @assert N > 2

p = zero(T)
c = v[1]
n = v[2]
@inbounds u[1] = d[1] * c + du[1] * n
@avx for i in 2:(N-1)
p = c
c = n
n = v[i+1]
u[i] = dl[i-1] * p + d[i] * c + du[i] * n
end
p = c
c = n
@inbounds u[N] = dl[N-1] * p + d[N] * c
u
end

let N = 5
T = Float64
A = Tridiagonal(rand(T, N-1), rand(T, N), rand(T, N-1))
v = rand(T, N)
u1 = Array{T}(undef, N)
u2 = Array{T}(undef, N)

matmul_avx!(u1, A, v)
mul!(u2, A, v)
u1 - u2
end

#+RESULTS:
: 5-element Vector{Float64}:
: 0.0
: 5.551115123125783e-17
: 0.3137201723219849
: 0.40281562647422675
: 0.29883752908665995
```
______

One thing I tried in order to fix this was to use an intermediate array to store `p, c, n`, but `@avx` didn't like that, giving an `undefvarerror`:
```julia
function matmul2_avx!(u::AbstractVector{T}, A::Tridiagonal{T}, v::AbstractVector{T}) where {T}
@assert length(u) == size(A,1) == size(A,2) == length(v)
dl, d, du = A.dl, A.d, A.du
N = length(u)
if N == 0
nothing
elseif N == 1
u[1] = d[1] * v[1]
return u
else
dep = MVector((zero(T), v[1], v[2]))
@inbounds u[1] = d[1] * dep[2] + du[1] * dep[3]
@avx for i in 2:(N-1)
dep[1] = dep[2]
dep[2] = dep[3]
dep[3] = v[i+1]
u[i] = dl[i-1] * dep[1] + d[i] * dep[2] + du[i] * dep[3]
end
dep[1] = dep[2]
dep[2] = dep[3]
@inbounds u[N] = dl[N-1] * dep[1] + d[N] * dep[2]
end
u
end

let N = 5
T = Float64
A = Tridiagonal(rand(T, N-1), rand(T, N), rand(T, N-1))
v = rand(T, N)
u1 = Array{T}(undef, N)
u2 = Array{T}(undef, N)

matmul2_avx!(u1, A, v)
mul!(u2, A, v)
u1 - u2
end

#+RESULTS:
:RESULTS:
# [goto error]
#+begin_example
UndefVarError: ######RHS###4######5### not defined

Stacktrace:
[1] matmul2_avx!(u::Vector{Float64}, A::Tridiagonal{Float64, Vector{Float64}}, v::Vector{Float64})
@ Main ./In[179]:41
[2] top-level scope
@ In[184]:8
[3] eval
@ ./boot.jl:360 [inlined]
[4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1094
#+end_example
:END:
```

Beitragsleitfaden

Für dieses Repository ist kein Beitragsleitfaden indexiert

Rechercherichtung

Start with the supplied matmul_avx! reproduction and the @avx loop, then compare it with the working matmul! version and the matmul2_avx! MVector attempt. Run the N=5 comparison against mul! and investigate the incorrect values and UndefVarError; done means the vectorized loop returns the same results without the error.

Vom Indexierungsmodell aus dem Issue-Text verfasst.

Bewertung

Tech-Stack
julia
Bereich
performance
Issue-Typ
Bug
Schwierigkeit
4/5
Geschätzter Aufwand
3-5 Tage
Aktivitätsstatus
Veraltet
Klarheit
Größtenteils klar
Anfängerfreundlichkeit
42/100

Neue Issues direkt in Ihr Postfach

Eine kurze Übersicht über anfängerfreundliche GitHub-Issues.