JuliaSIMD / JuliaSIMD/LoopVectorization.jl

`@avx` incorrectly handling variable redefinitions between iterations

Aperta
#266 5 commenti 0 reazioni 0 assegnatari Vedi su GitHub
Lingua principale
Julia
Stelle
789
Fork
73
Metriche di merge delle PR
Nessuna PR unita negli ultimi 30g

Descrizione

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:
```

Guida per i contributori

Nessuna guida per i contributori indicizzata per questo repository

Valutazione

Questa issue non è ancora stata valutata.

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.