[matthias314]

Since we now have fast iteration over FastContiguousSubArray (#48720), it could make sense to implement a special method for mapreduce for "reversed" FastSubArray, that is, for those step with value -1 (w2 below). Since for mapreduce we are allowed to process the array elements in any order, we can reverse the view to get something that has a much better chance of SIMD compiler optimization. Here are some benchmarks:

g(a) = reduce(+, a; init = zero(eltype(a)))
v = rand(Int8, 2^16);
w1 = @view v[begin:end];
w2 = @view v[end:-1:begin];

julia> @btime g($v); @btime g($w1); @btime g($w2)

# master
  574.189 ns (0 allocations: 0 bytes)
  641.550 ns (0 allocations: 0 bytes)
  834.255 μs (0 allocations: 0 bytes)
# this PR
  572.086 ns (0 allocations: 0 bytes)   # not affected by the PR
  675.769 ns (1 allocation: 48 bytes)   # around 30ns slower
  696.621 ns (1 allocation: 48 bytes)   # over 1000 times faster

The not implemented part of #48720 gave a speed-up factor of around 17 for FastSubArray{Int8}. Even if that is implemented later, this special method would still be around 70 times faster than that. Maybe this justifies to have a special method for this case. Maybe it doesn't.

Note: There is currently one issue with the PR. The new method allocates 48 bytes, which seems to take around 30ns (the time difference for w1 above). I don't know how to get rid of this, but I'm hoping that somebody else has an idea. I have the feeling that it is somehow related to the treatment of keyword arguments.

Also note that I'm using invoke to avoid any interference with the implementation of mapreduce for AbstractArray.

[N5N3]

This is unacceptable for op without commutability. (e.g. * for String, * for Matrix, vcat for Vector)

[matthias314]

Good point! In principle this is easy to fix -- just replace op by (x,y) -> op(y,x). Unfortunately, this seems to prevent the compiler from optimizing.

[mcabbott]

Or restrict to known commutative operations. IMO this should probably be hooked on later, to dispatch on dims rather than testing if T == Colon. Something similar is done for sums of Adjoint arrays, around here:

julia/stdlib/LinearAlgebra/src/adjtrans.jl

Lines 384 to 385 in 8068e44

	Base._mapreduce_dim(f, op::CommutativeOps, init::Base._InitialValue, A::Adjoint, dims::Colon) =
	Base._mapreduce_dim(f∘adjoint, op, init, parent(A), dims)

[matthias314]

Or restrict to known commutative operations.

I agree. The whole point of this PR is to enable SIMD optimization for views with step value -1. For SIMD, I can only think of BitInteger and float types and commutative operations (+, *, &, |, xor, maybe min and max) plus the non-commutative subtraction and division.

[N5N3]

just replace op by (x,y) -> op(y,x).

We have Base.FlipArgs

Unfortunately, this seems to prevent the compiler from optimizing.

This is a seperate problem. Once user provide init, Base.reduce always use foldl branch. (Thus disable many optimization)

[mikmoore]

Or restrict to known commutative operations.

To echo some above posts: I really dislike the idea of restricting to some specific set of types and operators. Such a special case has no support or extensibility for new operations. It's not just single primitive operations on primitive types that can SIMD, it's also compositions of those operations. Likewise, it can support other types that wrap these primitive types for slightly different behavior (which may or may not change the behavior on the suggested operations). It's also vulnerable to incorrect results when evil users do evil piracy, even if they've likely doomed themselves in other ways.

I'd much prefer we use a Base.FlipArgs(op) composition and let the compiler do its thing. If "the compiler doing its thing" is dumb, or we fall to foldl because of incomplete special-casing, that's an orthogonal issue.

[matthias314]

Here is a new version. I got it working for possibly non-commutative operations. As far as I can tell, it was the implementation of mapfoldr that prevented the compiler from optimizing. I've added a (naive?) mapfoldr method for AbstractArray. Please check that I didn't miss any case to consider there. Also, I've slightly rearranged the code for _mapreduce_dim to reduce the number of methods needed to avoid method ambiguities.

Here are some benchmarks: For a given type T, I define

v = rand(T, 2^16); w = @view(v[end:-1:begin]); i = zero(T)

Then I get

T = Int8	master	this PR
foldl(+,w)	817.218 μs	51.930 μs
foldr(+,w)	812.553 μs	666.411 ns
reduce(+,w)	22.043 μs	1.509 μs
reduce(+,w; init=i)	818.110 μs	669.516 ns

T = Int64	master	this PR
foldl(+,w)	823.908 μs	39.416 μs
foldr(+,w)	833.575 μs	8.162 μs
reduce(+,w)	23.128 μs	8.576 μs
reduce(+,w; init=i)	813.772 μs	7.941 μs

T = Float64	master	this PR
foldl(+,w)	830.546 μs	69.227 μs
foldr(+,w)	835.472 μs	7.729 μs
reduce(+,w)	67.846 μs	8.586 μs
reduce(+,w; init=i)	839.744 μs	7.947 μs

For the last table I used @fastmath. I tested all BitInteger types with +, *, &, |, xor, min, max and Float16, Float32 and Float64 with +, *, min, max. The results are similar. I didn't see any allocations nor any significant change in runtime compared to master when using the original vector v or view(v, :) instead of w. Well, in a few cases I got differences for v and foldr by a factor of up to 2, but I couldn't reproduce that.

[matthias314]

@N5N3 Do you think this is good now? In that case I would rebase the commits into a single one.

A related issue (outside the scope of this PR, I think) would be to understand why reduce without init is so much faster than foldl and foldr. Getting that speed for folding would dramatically improve performance for FastSubArray with stride1 different from 1 and -1.

EDIT: This is because _mapreduce uses @inbounds. If that's allowed, then we could speed up things easily!

[N5N3]

No need to change these?

[matthias314]

Without this change, I would have to define two separate methods for _mapreduce_dim for FastSubArray, one with _InitialValue and one without. Moreover, the method with _InitialValue would have to go to reducedim.jl instead of subarray.jl because _InitialValue is not yet defined when subarray.jl is read. I wanted to keep all the SubArray code in subarray.jl.

However, with your idea to change _mapreduce instead, no change would be necessary in reducedim.jl.

[N5N3]

Suggested change

	function mapfoldl_impl(f::F, op::OP, init::I, a::FastSubArray) where {F,OP,I}
	if a.stride1 == -1
	b = reverse_fastsubarray_1(a)
	@invoke mapfoldr_impl(f::F, FlipArgs(op)::FlipArgs{OP}, init::I, b::AbstractArray)
	else
	@invoke mapfoldl_impl(f::F, op::OP, init::I, a::AbstractArray)
	end
	end
	function mapfoldr_impl(f::F, op::OP, init::I, a::FastSubArray) where {F,OP,I}
	if a.stride1 == -1
	b = reverse_fastsubarray_1(a)
	@invoke mapfoldl_impl(f::F, FlipArgs(op)::FlipArgs{OP}, init::I, b::AbstractArray)
	else
	@invoke mapfoldr_impl(f::F, op::OP, init::I, a::AbstractArray)
	end
	end
	function mapfoldl_impl(f::F, op::OP, nt, a::FastSubArray) where {F,OP}
	op′, itr′ = _xfadjoint(BottomRF(op), Generator(f, a))
	if itr′ isa FastSubArray && itr′.stride1 == -1
	return foldl_impl(op′, nt, _reverse_iter(reverse_fastsubarray_1(itr′)))
	else
	return foldl_impl(op′, nt, itr′)
	end
	end
	function mapfoldr_impl(f::F, op::OP, nt, a::FastSubArray) where {F,OP}
	op′, itr′ = _xfadjoint(BottomRF(FlipArgs(op)), Generator(f, a))
	if itr′ isa FastSubArray && itr′.stride1 == -1
	return foldl_impl(op′, nt, reverse_fastsubarray_1(itr′))
	else
	return foldl_impl(op′, nt, _reverse_iter(itr′))
	end
	end

Perhaps this is better.

[matthias314]

Your proposal copies the logic from mapfoldl_impl and mapfoldr_impl. I didn't want to interfere with how these functions are actually implement for AbstractArray because any future change of that would have to be mirrored for FastSubArray. Does this convince you? Otherwise I change as you have suggested.

[N5N3]

Does this convince you?

Make sense. I just think a FlipArgs(FlipArgs(...)) is ugly and the current implement is wrong as FlipArgs(op) !<: FlipArgs{OP} if op isa constructor, e,g, Pair.

In fact, when I switch to the above code, I increasingly feel that this is not a good idea.
If we just want to avoid boundcheck, then I think an overall fix for #48773 would be much better.
(If guys still dislike @propagte_inbounds, then we can add a _foldl_impl for AbstractArray.)

If we want more auto-vectorlization, I don't think stride1 == -1 is that special. And we still need tools for other simdable stride1.

[matthias314]

FlipArgs(op) !<: FlipArgs{OP} if op is a constructor, e,g, Pair

We absolutely have to avoid wrong code. However, I get

julia> using Base: FlipArgs
julia> op = Pair; FlipArgs(op) isa FlipArgs{typeof(op)}
true
julia> op = FlipArgs(+); FlipArgs(op) isa FlipArgs{typeof(op)}
true

an overall fix for #48773 would be much better

I fully agree. I'd like to see what comes out of my Discourse post. So far there has been no clear statement whether @inbounds is acceptable inside folding operations.

If we want more auto-vectorization, I don't think stride1 == -1 is that special. And we still need tools for other simdable stride1.

Is there auto-vectorization for non-contiguous data? That would be great. However, I would expect that even in that case the step value has to known at compile time. If that's true, then there is a limit to what we can do. I'm thinking that the -1 case is still fairly common. But again, a general solution would be much better.

[N5N3]

We absolutely have to avoid wrong code. However, I get

My fault, didn't notice that FlipArgs(op) is not stable for constructor. Perhaps we should fix it.

I would expect that even in that case the step value has to known at compile time.

You are right, without compile time info LLVM would generate a simd branch for dense input and check if it's valid at runtime.
Since this only benifits simdable cases, it makes sense to me to to ask users to provide more static information if this does matter.
We already support range with static step in Static.jl, but that's not enough as we use the cached stride1 in getindex.

[N5N3]

I think we'd better overload _mapreduce(f, op, ::IndexLinear, ::FastSubArray) instead.

[matthias314]

That would slow mapreduce with given init down because _mapreduce only covers the case of mapreduce without init. If init is given, then mapreduce calls mapfoldl_impl. For a "reversed" view, mapfoldl_impl is translated into mapfoldr_impl of the parent array. The current PR instead uses mapfoldl_impl of the parent array, which is 5 to 75 times faster, see the difference between foldl and foldr in the benchmarks above.

[N5N3]

But we have add mapfoldl(r)_impl specializaiton for FastSubArray above?

[matthias314]

Yes, that's exactly the point: For w = view(v, 3:-1:1) the following happens with your suggestion:

mapreduce(f, op, w; init = i) =
mapfoldl_impl(f, op, i, w) =
mapfoldr_impl(f, FlipArgs(op), i, view(v, 1:3))

For such a view (with known step value 1) mapfoldr is much slower than mapfoldl. It's the same as for the underlying array:

julia> v = rand(Int8, 2^16); w = @view v[begin:end];
julia> @btime foldl(+, $v);
  578.033 ns (0 allocations: 0 bytes)
julia> @btime foldl(+, $w);
  644.651 ns (0 allocations: 0 bytes)
julia> @btime foldr(+, $v);
  34.613 μs (0 allocations: 0 bytes)
julia> @btime foldr(+, $w);
  39.116 μs (0 allocations: 0 bytes)

EDIT: For comparison, I've created the branch matthias314/julia@ea64de2 where @inbounds is used for mapfoldl and friends and also in the lowering of for loops. There one gets:

julia> @btime foldl(+, $v);
  574.355 ns (0 allocations: 0 bytes)
julia> @btime foldl(+, $w);
  574.304 ns (0 allocations: 0 bytes)
julia> @btime foldr(+, $v);
  1.114 μs (0 allocations: 0 bytes)
julia> @btime foldr(+, $w);
  1.108 μs (0 allocations: 0 bytes)

[DilumAluthge]

This PR has had no activity for two years+, so I'll close it as stale.

@matthias314 If you're still interested in working on this PR, let me know, and I can re-open it.