SIMD median

SIMDmedian.jl

using SIMD

Base.@pure simdwidth(::Type{T}) where {T} = Int(256/8/sizeof(T))

@inline function median3(a,b,c)
    max(min(a,b), min(c,max(a,b)))
end
@inline function median5(a,b,c,d,e)
    # https://stackoverflow.com/questions/480960/code-to-calculate-median-of-five-in-c-sharp
    f=max(min(a,b),min(c,d))
    g=min(max(a,b),max(c,d))
    median3(e,f,g)
end

@noinline function median5_vectors!(out, a,b,c,d,e)
    K = simdwidth(eltype(out))
    N = length(out)
    T = eltype(out)
    V = Vec{K,T}
    @assert mod(N,K) == 0

    @inbounds for i in 1:K:N
        va = vload(V,a, i)
        vb = vload(V,b, i)
        vc = vload(V,c, i)
        vd = vload(V,d, i)
        ve = vload(V,e, i)
        vo = median5(va,vb,vc,vd,ve)
        vstore(vo,out, i)
    end
    out
end

using BenchmarkTools
T = UInt8
N = 10^6
N = N ÷ simdwidth(T) * simdwidth(T)
out, a,b,c,d,e = [rand(T,N) for _ in 1:6]
@benchmark median5_vectors!(out, a,b,c,d,e)