Extrema is slower than maximum + minimum

[mbauman]

On AVX-512 systems, calling extrema(x) is slower than calling maximum and then minimum. Looks like we need to incorporate the same techniques from #30320 over to extrema in order to get it to vectorize in a similar manner.

@mbauman Can you provide the test case you used? I ran the following code:

arrs = []

for i in 1:100000
    insert!(arrs, 1, rand(1000))
end

println("running extrema")
@time map(arrs) do arr
    extrema(arr)
end

println("running minimum, maximum")
@time map(arrs) do arr
    (minimum(arr), maximum(arr))
end

and got this output:

running extrema
  0.520760 seconds (301.98 k allocations: 14.811 MiB)
running minimum, maximum
  0.741328 seconds (279.30 k allocations: 13.814 MiB)

julia version:

julia version 1.1.0

[mbauman]

Yes, that test case reproduces the behavior on master for me. It won't on 1.1.0 as the PR that improved maximum and minimum was after 1.1.0. Of course, note that you'll also need to have a relatively recent computer to see the difference.

[nlw0]

I got a slower extrema in my Core i5 machine.

running extrema
  314.316 ms (100003 allocations: 4.58 MiB)
running minimum, maximum
  224.977 ms (100003 allocations: 4.58 MiB)

[nlw0]

I would like to try working on this, and it would be great to hear what kind of patch is expected. I was thinking maybe mapreduce_impl could be turned into a generated function that can handle multiple functions simultaneously. Does that sound good?

One question I have is about the unrolling. Ideally this should be adapted for each architecture and data type. How much of this should be up to LLVM or to Julia? Should we prioritize generating code that the LLVM can vectorize? The alternative I imagine would be controlling parameters like this unrolling depending on the architecture.

[nlw0]

Thinking a bit more about this, maybe the best of worlds would be if this function could simply be implemented with a general fold-left, and it should result in a good vectorized machine code relying on the LLVM automatic vectorization. Should we concentrate our efforts on making sure that a generic fold works better in more cases, or should we first make sure that a few specific functions have high-performance implementations?

[mbauman]

I was thinking of the much simpler option of copy-pasting the techniques used in maximum/minimum (here:

julia/base/reduce.jl

Lines 481 to 518 in 2b49d26

	a1 = @inbounds A[first]
	v1 = mapreduce_first(f, op, a1)
	v2 = v3 = v4 = v1
	chunk_len = 256
	start = first + 1
	simdstop = start + chunk_len - 4
	while simdstop <= last - 3
	# short circuit in case of NaN
	v1 == v1 || return v1
	v2 == v2 || return v2
	v3 == v3 || return v3
	v4 == v4 || return v4
	@inbounds for i in start:4:simdstop
	v1 = _fast(op, v1, f(A[i+0]))
	v2 = _fast(op, v2, f(A[i+1]))
	v3 = _fast(op, v3, f(A[i+2]))
	v4 = _fast(op, v4, f(A[i+3]))
	end
	checkbounds(A, simdstop+3)
	start += chunk_len
	simdstop += chunk_len
	end
	v = op(op(v1,v2),op(v3,v4))
	for i in start:last
	@inbounds ai = A[i]
	v = op(v, f(ai))
	end
	# enforce correct order of 0.0 and -0.0
	# e.g. maximum([0.0, -0.0]) === 0.0
	# should hold
	if isbadzero(op, v)
	for i in first:last
	x = @inbounds A[i]
	isgoodzero(op,x) && return x
	end
	end
	return v

) into the body of the extrema implementation (here:

julia/base/operators.jl

Lines 471 to 483 in 2b49d26

	y = iterate(itr)
	y === nothing && throw(ArgumentError("collection must be non-empty"))
	(v, s) = y
	vmin = vmax = f(v)
	while true
	y = iterate(itr, s)
	y === nothing && break
	(x, s) = y
	fx = f(x)
	vmax = max(fx, vmax)
	vmin = min(fx, vmin)
	end
	return (vmin, vmax)

).

Revamping mapreduce to support multiple functions is a much bigger task (and I'm not sure we want to add that complication to the already complex mapreduce machinery, either), and while improving LLVM's automatic vectorization is always a worthy goal it's particularly challenging in cases like these.

[nlw0]

I am still trying to understand why the specialization of mapreduce_impl was needed in the first place, and what are the symmetries between minimum and extrema. I am reading the history of the file and I cannot make complete sense of it.

The generic mapreduce_impl contains a @simd for line, and produces vectorized code to op=max from what I can tell. The only reason I see for the specialized function is the short-circuiting. Are we interested in having this on extrema too? Also, these specialized functions are not using @simd for, but a while instead. Is this on purpose?

This idea I had at first of supporting multiple functions doesn't really make sense, sorry. I am thinking now on something much simpler, which would be to implement this function as a fold like:
myextrema(arr) = Base.mapfoldl(identity, ((mi,ma),v)->(min(mi,v),max(ma,v)), arr[2:end], init=(arr[1],arr[1]))

Which of course we can specialize if it makes sense. What I see now is that this doesn't seem to vectorize, and I wonder if just refactoring mapfoldl_impl would do it, and that is what I see as something that might be nice to investigate, since it might also improve other functions based on it.

It doesn't seem to matter for this specific function, though, because from what I can tell the current extrema is already vectorized. So I suppose in the end this is just about using _min_fast and _max_fast?

[nlw0]

I quickly tried to use the _min_fast functions and saw no change. And the machine code is already showing the vcmpordsd instruction. Maybe this will really require something like that unrolling?

[mbauman]

Note that just seeing vector instructions in @code_native set doesn't actually mean the processor is using the parallelism of the SIMD unit... it can also just be using SIMD instructions with single values. For example, take a look at @code_native max(1.0, 2.0) — that'll show a whole slew of "SIMD" instructions (including that vcmpordsd) operating on that single data. It's much easier to look at @code_llvm to see if things are vectorizing.

For example, you can see that extrema actually does vectorize nicely for integers:

julia> @code_llvm extrema(rand(Int, 1000))
# ... contains things like:
vector.ph:                                        ; preds = %L40.lr.ph.L40.lr.ph.split_crit_edge
      %n.vec = and i64 %27, -32
      %ind.end = or i64 %n.vec, 1
      %ind.end26 = or i64 %n.vec, 2
      %minmax.ident.splatinsert = insertelement <8 x i64> undef, i64 %24, i32 0
      %minmax.ident.splat = shufflevector <8 x i64> %minmax.ident.splatinsert, <8 x i64> undef, <8 x i32> zeroinitializer
      %minmax.ident.splatinsert59 = insertelement <8 x i64> undef, i64 %24, i32 0
      %minmax.ident.splat60 = shufflevector <8 x i64> %minmax.ident.splatinsert59, <8 x i64> undef, <8 x i32> zeroinitializer
      br label %vector.body

vector.body:                                      ; preds = %vector.body, %vector.ph
      %index = phi i64 [ 0, %vector.ph ], [ %index.next, %vector.body ]
      %vec.phi = phi <8 x i64> [ %minmax.ident.splat, %vector.ph ], [ %48, %vector.body ]
# ...

Key things to see are the vector.body label, all the <8 x i64> vectors. If you look at @code_llvm extrema(rand(1000)) with floating point values, though, you'll see no vector.body, no <N x double> vectors, and lots of plain old double comparisons.

[chriselrod]

vcmpordsd isn't a SIMD instruction. sd means "single double". If it was named pd for "packed double", it would be SIMD.