JuliaSIMD / JuliaSIMD/LoopVectorization.jl
`@avx` incorrectly handling variable redefinitions between iterations
- 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