On master, mapreduce calls mapfoldl for inputs with IndexCartesian style, which brings performance pain if the 1st dimension is long enough for vectorlization.

This PR implements pairwise mapreduce_impl for better performance by cutting the highest splitable dimension in half.
(As we can't reorder the inputs, and CartesianPartition has much higher overhead)

Test added.

A simple benchmark

julia> a = randn(100, 100);

julia> slow(a) = view(a, axes(a)...)
slow (generic function with 1 method)

julia> @btime reduce(+, slow($a))
  1.870 μs (0 allocations: 0 bytes)
88.67608098492965

julia> @btime foldl(+, slow($a))
  9.700 μs (0 allocations: 0 bytes)
88.67608098492977