JuliaSIMD / JuliaSIMD/LoopVectorization.jl

Performance of convolutions over 3-tuples

Abierto

#342 abierto el 25/9/2021 5 comentarios 0 reacciones 0 asignados Ver en GitHub

Lenguaje dominante: Julia
Estrellas: 789
Forks: 73
Métricas de merge de PR: Sin PR fusionados en 30 d

Descripción

Ref: https://github.com/JuliaSIMD/LoopVectorization.jl/pull/341

@timholy (Feel free to also share your benchmark code)

I can run benchmarks on a few different computers.
I also own Skylake-X/Cascadelake (Skylake-era AVX512 with 2x FMA), Haswell (AVX2), and Apple M1.

For now, I am running on:
```julia
julia> versioninfo()
Julia Version 1.8.0-DEV.593
Commit 53aa65ad3c* (2021-09-24 19:45 UTC)
Platform Info:
OS: Linux (x86_64-redhat-linux)
CPU: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
```
Which has AVX512 with 1x FMA on port 0.
In addition, it also has vector shuffles on port 5.

Contrast this with my Skylake-X/Cascadelake CPUs, which have vector shuffles on port 5, and FMAs on both ports 0 and 5. What this means is that on Skylake-X/Cascadelake, performing vector shuffles is using resources that could otherwise be spent on FMA.
But on Tiger Lake, they're independent / fully parallel.

So, this means on Tiger Lake, we would ideally see comparable flops between greyscale (requires vector moves and fma) and array-of-structs [colors] (requires vector moves, shuffles, and fma).

Functions:
```julia
using LoopVectorization, OffsetArrays, BenchmarkTools
function smoothdim_kernel_tile!(out, z, src::AbstractArray, kernel::AbstractVector, Rpre::CartesianIndices, axout_tile, Rpost::CartesianIndices, axkernel = axes(kernel, 1))
# z is an "overflow-safe zero". This is needed in cases where both `src` and `kernel` have eltypes vulnerable to arithmetic overflow.
@fastmath @inbounds for Ipost in Rpost
for i in axout_tile
for Ipre in Rpre
tmp = convert(eltype(out), z)
for j in axkernel
tmp += oftype(z, src[Tuple(Ipre)...,i+j,Tuple(Ipost)...])*kernel[j]
end
out[Tuple(Ipre)...,i,Tuple(Ipost)...] = tmp
end
end
end
out
end
function smoothdim_kernel_tile_sa!(out, z, src::AbstractArray, kernel::AbstractVector, Rpre::CartesianIndices, axout_tile, Rpost::CartesianIndices, axkernel = LoopVectorization.ArrayInterface.static(-6):LoopVectorization.ArrayInterface.static(6))
# z is an "overflow-safe zero". This is needed in cases where both `src` and `kernel` have eltypes vulnerable to arithmetic overflow.
@fastmath @inbounds for Ipost in Rpost
for i in axout_tile
for Ipre in Rpre
tmp = convert(eltype(out), z)
@simd ivdep for j in axkernel
tmp += oftype(z, src[Tuple(Ipre)...,i+j,Tuple(Ipost)...])*kernel[j]
end
out[Tuple(Ipre)...,i,Tuple(Ipost)...] = tmp
end
end
end
out
end
function smoothdim_kernel_tile_avx_inner!(out, z, src::AbstractArray, kernel::AbstractVector, Rpre::CartesianIndices, axout_tile, Rpost::CartesianIndices)
axkernel = axes(kernel, 1)
zout = convert(eltype(out), z)
for Ipost in Rpost
for i in axout_tile
@turbo for Ipre in Rpre
tmp = zout
# tmp = convert(eltype(out), z) # failing to hoist this leads to an "UndefVarError: tmp not defined"
for j in axkernel
tmp += oftype(z, src[Ipre,i+j,Ipost])*kernel[j]
end
out[Ipre,i,Ipost] = tmp
end
end
end
out
end
function smoothdim_kernel_tile_avx_middle!(out, z, src::AbstractArray, kernel::AbstractVector, Rpre::CartesianIndices, axout_tile, Rpost::CartesianIndices)
axkernel = axes(kernel, 1)
zout = convert(eltype(out), z)
for Ipost in Rpost
@turbo for i in axout_tile
for Ipre in Rpre
tmp = zout
# tmp = convert(eltype(out), z) # failing to hoist this leads to an "UndefVarError: tmp not defined"
for j in axkernel
tmp += oftype(z, src[Ipre,i+j,Ipost])*kernel[j]
end
out[Ipre,i,Ipost] = tmp
end
end
end
out
end
function smoothdim_kernel_tile_avx_outer!(out, z, src::AbstractArray, kernel::AbstractVector, Rpre::CartesianIndices, axout_tile, Rpost::CartesianIndices)
axkernel = axes(kernel, 1)
zout = convert(eltype(out), z)
@turbo for Ipost in Rpost
for i in axout_tile
for Ipre in Rpre
tmp = zout
# tmp = convert(eltype(out), z) # failing to hoist this leads to an "UndefVarError: tmp not defined"
for j in axkernel
tmp += oftype(z, src[Ipre,i+j,Ipost])*kernel[j]
end
out[Ipre,i,Ipost] = tmp
end
end
end
out
end
```
The reason for the `Tuple(Ipost)...` type indexing is that this was needed for type stability with 5-d arrays.
Note that the middle and inner variants of the `@turbo` functions need some work for type stability, too, but this should be fixable within LoopVectorization.

This is probably a new regression in Julia 1.8.

Benchmarks code for `d=1`, `Float32` and then `NTuple{3,Float32}`:
```julia
T = Float32; M = 11; σ = 3;
ax = OffsetArray(-6:6, -6:6);
kernel = float(T)[exp(-i^2/(2*σ^2)) for i in ax]; # 1-d Gaussian kernel
kernel ./= sum(kernel);
srcax = 1+firstindex(ax):M+lastindex(ax);

Mpad = length(srcax);
src = OffsetArray(rand(T,Mpad,Mpad,Mpad),srcax,srcax,srcax);
dest1 = Array{float(T),3}(undef,M,M,M);
dest2 = similar(dest1);

d=1;
Rpre = CartesianIndices(axes(dest1)[1:d-1]);
Rpost = CartesianIndices(axes(dest1)[d+1:end]);
smoothdim_kernel_tile!( dest1, float(zero(T)), src, kernel, Rpre, axes(dest1, d), Rpost);
smoothdim_kernel_tile_avx_inner!(dest2, float(zero(T)), src, kernel, Rpre, axes(dest2, d), Rpost);
dest1 ≈ dest2

@benchmark smoothdim_kernel_tile!( $dest1, $(float(zero(T))), $src, $kernel, $Rpre, axes($dest1, $d), $Rpost)
@benchmark smoothdim_kernel_tile_sa!( $dest1, $(float(zero(T))), $src, $kernel, $Rpre, axes($dest1, $d), $Rpost)
@benchmark smoothdim_kernel_tile_avx_inner!($dest2, $(float(zero(T))), $src, $kernel, $Rpre, axes($dest2, $d), $Rpost)
@benchmark smoothdim_kernel_tile_avx_middle!($dest2, $(float(zero(T))), $src, $kernel, $Rpre, axes($dest2, $d), $Rpost)
@benchmark smoothdim_kernel_tile_avx_outer!($dest2, $(float(zero(T))), $src, $kernel, $Rpre, axes($dest2, $d), $Rpost)

src_t = OffsetArray([(rand(T),rand(T),rand(T)) for i in srcax, j in srcax, k in srcax], srcax, srcax, srcax);
srcr = reinterpret(reshape, T, src_t);
fT = float(T);
dest1_t = Array{Tuple{fT,fT,fT}}(undef, M, M, M);
dest2_t = similar(dest1_t);
dest1r = reinterpret(reshape, fT, dest1_t);
dest2r = reinterpret(reshape, fT, dest2_t);

d = 1;
Rprer = CartesianIndices((LoopVectorization.Static(1):LoopVectorization.Static(3), axes(dest1_t)[1:d-1]...));
Rpostr = CartesianIndices(axes(dest1_t)[d+1:end]);
smoothdim_kernel_tile!( dest1r, float(zero(T)), srcr, kernel, Rprer, axes(dest1_t, d), Rpostr);
smoothdim_kernel_tile_avx_inner!(dest2r, float(zero(T)), srcr, kernel, Rprer, axes(dest1_t, d), Rpostr);

@benchmark smoothdim_kernel_tile!( $dest1r, float(zero($T)), $srcr, $kernel, $Rprer, axes($dest1_t, $d), $Rpostr)
@benchmark smoothdim_kernel_tile_sa!( $dest1r, float(zero($T)), $srcr, $kernel, $Rprer, axes($dest1_t, $d), $Rpostr)
@benchmark smoothdim_kernel_tile_avx_inner!($dest2r, float(zero($T)), $srcr, $kernel, $Rprer, axes($dest2_t, $d), $Rpostr)
@benchmark smoothdim_kernel_tile_avx_middle!($dest2r, float(zero($T)), $srcr, $kernel, $Rprer, axes($dest2_t, $d), $Rpostr)
@benchmark smoothdim_kernel_tile_avx_outer!($dest2r, float(zero($T)), $srcr, $kernel, $Rprer, axes($dest2_t, $d), $Rpostr)
```

The `sa` variant is using a static axis and `@simd ivdep` to encourage LLVM to SIMD.

Among the `@turbo` variants, the `inner` are like those from the tests.
The `middle` and `outer` move the `@turbo` to outer loops, giving LV more options for optimization.
This is particularly relevant for `d=1`, where most of the loops are outer-loops, meaning LV doesn't have many options when `@turbo` is applied toward the inside.
In fact, when `d=1`, `Rpre` is 0-dimensional, meaning the only loop LV can optimize is the inner-most loop.

Benchmark results for `Float32`:
```julia
julia> @benchmark smoothdim_kernel_tile!( $dest1, $(float(zero(T))), $src, $kernel, $Rpre, axes($dest1, $d), $Rpost)
BenchmarkTools.Trial: 10000 samples with 6 evaluations.
Range (min … max): 5.661 μs … 10.469 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 5.753 μs ┊ GC (median): 0.00%
Time (mean ± σ): 5.888 μs ± 284.296 ns ┊ GC (mean ± σ): 0.00% ± 0.00%

▆▆ ▃█▆ ▃▆▁ ▂▄▁ ▂▃ ▂▇▄ ▂
███▁███▅███▆▅███▇███▆▇▇▇▆▇▆▅▅▆▆▆▆▃▅▄▆▅▅▅▃▇▆▅▅███▇▆▆▇▆▇▄▅▅▇▇ █
5.66 μs Histogram: log(frequency) by time 6.56 μs <