Core: Use four accumulators in the vectorized reductions

redux_impl<Func, Evaluator, LinearVectorizedTraversal, NoUnrolling> carries two packet accumulators. That makes the loop run at the latency of func.packetOp rather than at its throughput: with a four-cycle vaddpd/vmaxpd, two accumulators cap the loop at four cycles per two packets no matter how many execution ports are idle.

This MR uses four independent accumulators instead. The unrolled body is bounded by numext::round_down(size - alignedStart, 4 * packetSize), matching the packetSize bound on the line below it, and the four accumulators merge as (res0 + res1) + (res2 + res3), two packetOp latencies deep rather than three.

Short lengths keep the cheap shapes they had. Two and three packets take a straight-line path with none of the trip-count setup a loop needs, and the one to three packets the unrolled body leaves over are folded into the three accumulators that are still independent, before the merge rather than after it. The first revision of this MR gave all of that to one shared tail loop and lost up to 1.55x on two- to seven-packet reductions; the second commit is that fix, measured below.

The change alters a dependency chain, not code generation, so it needs no compiler or ISA gate.

redux_impl<Func, Evaluator, SliceVectorizedTraversal, Unrolling> had the same problem in a worse form: it feeds one packet_res through the whole (j, i) traversal, so m.block(...).sum() or m.topLeftCorner(...).maxCoeff() runs at one packetOp latency per packet. It gets four accumulators too, merged the same way, and where they come from depends on the panel:

• four packets or more per panel: each accumulator takes one packet of every group of four, and the one to three packets the panel has left over go into the three that are still independent;
• one to three packets per panel: a panel cannot supply four independent packets, so the independence comes from the outer dimension instead and each accumulator takes its own panel;
• fewer than four such panels: at most nine packets in total, which keeps the single accumulator.

All three live in run() alongside one copy of the predux and the scalar tail, the same shape the linear specialization has. An earlier revision put them in a packetRedux() helper; since nothing forced that inline, GCC sometimes left it as a call and spilled the accumulators across it, which cost up to 1.45x on narrow panels.

Measurements

benchmarks/Core/bench_reductions.cpp built against master and against this branch with identical flags -- the table below from the file as master has it, the two tables after it from the sizes and cases this MR adds to it. 13th Gen Intel(R) Core(TM) i7-13700HX, WSL2, taskset -c 2, GCC 13.3, -O3 -DNDEBUG -std=c++14. Each case was measured in its own process (see the note on methodology below); three master/branch pairs, seven repetitions each, median CPU ns.

Benchmark	ISA	master (ns)	this MR (ns)	Ratios (3 pairs)
VectorSum_float/4096	AVX2	115.41	66.67	1.73 1.71 1.75
VectorSum_double/4096	AVX2	232.02	128.02	1.84 1.81 1.79
VectorMaxCoeff_float/65536	AVX2	3910.11	1964.38	2.00 2.04 1.94
VectorMaxCoeff_double/65536	AVX2	8010.64	4149.83	1.88 1.97 1.94
VectorSquaredNorm_float/65536	AVX2	3995.11	1961.93	2.11 2.03 1.97
VectorSquaredNorm_double/16384	AVX2	1861.32	899.44	2.18 2.09 1.95
VectorNorm_float/16384	AVX2	885.00	466.88	1.92 1.81 1.96
VectorProd_double/4096	AVX2	420.55	220.84	1.97 1.83 1.91
VectorLpNormInf_double/16384	AVX2	2143.06	1115.29	1.68 1.92 2.21
MatrixNorm_double/256	AVX2	7304.44	3876.26	1.95 1.80 1.90
VectorMean_double/4096	AVX2	225.01	119.28	1.92 1.88 1.85
VectorMinCoeff_float/262144	AVX2	17433.42	9224.84	2.02 1.53 2.25
VectorSum_float/4096	SSE2	218.97	115.12	1.82 1.88 2.01
VectorSum_double/4096	SSE2	437.54	228.01	1.95 1.95 1.86
VectorMaxCoeff_float/65536	SSE2	7264.73	3693.19	1.98 1.93 1.99
VectorMaxCoeff_double/65536	SSE2	14698.15	7537.23	1.99 1.93 1.93
VectorSquaredNorm_float/65536	SSE2	3721.02	2633.68	1.43 1.45 1.36
VectorSquaredNorm_double/16384	SSE2	1841.86	1288.63	1.42 1.43 1.43
VectorNorm_float/16384	SSE2	910.75	639.91	1.44 1.40 1.43
VectorProd_double/4096	SSE2	857.79	446.49	1.88 1.98 1.91
VectorLpNormInf_double/16384	SSE2	4460.64	2246.26	1.96 1.94 2.06
MatrixNorm_double/256	SSE2	7494.44	5385.79	1.38 1.42 1.37
VectorMean_double/4096	SSE2	447.95	224.49	2.03 1.97 1.99
VectorMinCoeff_float/262144	SSE2	29776.19	15294.33	2.01 1.88 1.95

Short vectors

The sizes the second commit is about, in packets of the scalar. Same host and flags, BM_VectorSum with only that benchmark registered, each case in its own process, four alternating rounds of nine repetitions, --benchmark_min_time=0.02s. 4c9fc257 is the merge base, b8c5953a this MR's first commit, cd515583 its second. Median CPU ns:

Benchmark	packets	4c9fc257	b8c5953a	cd515583	first/base	second/base
VectorSum_double/8	2	0.745	1.052	0.734	1.41x	0.99x
VectorSum_double/12	3	0.756	1.171	0.741	1.55x	0.98x
VectorSum_double/24	6	1.071	1.456	1.035	1.36x	0.97x
VectorSum_double/28	7	1.166	1.475	1.049	1.27x	0.90x
VectorSum_float/16	2	0.815	0.972	0.797	1.19x	0.98x
VectorSum_float/24	3	0.910	1.270	0.857	1.40x	0.94x
VectorSum_float/48	6	1.237	1.524	1.180	1.23x	0.95x
VectorSum_float/56	7	1.338	1.504	1.318	1.12x	0.99x
VectorSum_double/4096	1024	206.960	110.963	107.789	0.54x	0.52x

-msse2 tells the same story at its own packet widths: 1.15-1.43x of the baseline after the first commit, 0.95-1.01x after the second, with VectorSum_double/16384 unchanged at 0.66-0.68x. Clang 18.1.3 already compiled the shared tail loop well enough to sit near the baseline, and gains 0.74-0.85x on the same short cases. Long reductions do not move: VectorSum_double/65536 1.02x, VectorSum_float/16384 1.00x, MatrixNorm_double/256 0.98x, all second commit over first.

bench_reductions.cpp could not see any of this, because VECTOR_SIZES started at 64 elements. It now takes the scalar as a parameter and registers two, three, six and seven packets of it ahead of the bandwidth-bound sizes, so every type and ISA brackets the same boundaries.

Block reductions

SliceVectorizedTraversal, which nothing in bench_reductions.cpp measured before this MR. BlockSum, BlockMaxCoeff and BlockSquaredNorm reduce a rows x cols block of a matrix eight larger in both directions, with the panel height given in packets of the scalar. Same host, flags and methodology; ratios are this MR over master:

Benchmark	panel	ratio	Benchmark	panel	ratio
BlockSum_double/4/8	1 packet	0.40x	BlockSum_float/8/8	1 packet	0.30x
BlockSum_double/8/64	2	0.48x	BlockSum_float/16/64	2	0.59x
BlockSum_double/12/64	3	0.43x	BlockSum_float/48/64	6	0.59x
BlockSum_double/16/64	4	0.74x	BlockSum_float/64/512	8	0.48x
BlockSum_double/24/64	6	0.56x	BlockSum_float/256/256	32	0.41x
BlockSum_double/32/512	8	0.48x	BlockMaxCoeff_float/16/64	2	0.35x
BlockSum_double/64/512	16	0.43x	BlockMaxCoeff_double/8/64	2	0.38x
BlockSum_double/256/256	64	0.40x	BlockMaxCoeff_double/64/512	16	0.26x
BlockSum_double/1024/1024	256	0.92x	BlockSquaredNorm_double/64/512	16	0.52x

The 1024x1024 case is DRAM-bound and keeps 1.1x of the 2.5x the cached sizes get.

Over the whole 438-case suite run end to end (two master/branch pairs, 7 repetitions, --benchmark_min_time=0.01s): 254 cases improved by more than 5% in both runs on AVX2 and 284 on SSE2.

Why the whole-suite run cannot be read case by case, and what the apparent regressions turn out to be

The same 438-case sweep reports 9 cases (AVX2) and 10 cases (SSE2) as more than 5% slower in both runs. Re-measuring each of them in its own process removes every one.

The AVX2 list is all working sets of 4-8 MB — VectorMean_float/1048576, VectorMinCoeff_double/1048576, MatrixNorm_double/1024 and so on. In isolation they land at 1.01-1.16 with per-pair ratios that swing from 0.32 to 1.26 on the same binary, i.e. they are DRAM-bandwidth-bound and dominated by host noise, not by this change.

Three SSE2 entries are more interesting because they reproduce: VectorReduxCommutativeOp_float/128 (0.72-0.74), VectorReduxUserOp_double/24 (0.69-0.70) and StridedRowSum_double/24 (0.75-0.80). None of them can execute the modified code. UserSumOp and CommutativeUserSumOp leave functor_traits::PacketAccess at its default false, so redux_traits never selects LinearVectorizedTraversal for them, and m.row(32) on a column-major matrix has a non-unit inner stride, so it has no packet path either.

Rebuilding just those three benchmarks into a three-case translation unit drops VectorReduxCommutativeOp_float/128 to 0.97-1.00. The other two persist at 0.82-0.89, so I disassembled them: BM_StridedRowSum<double> is 251 instructions in both binaries, differing only in branch target addresses and rip-relative displacements. The function is byte-identical and merely relocated. What is being measured on a 24-element scalar reduction that takes 3 ns is code placement — alignment and branch-predictor aliasing — which a one-line change anywhere in Redux.h will move either way.

This is the same effect that makes the end-to-end sweep unreliable per case, so I have reported the headline table from isolated runs only.

Validation

redux (all 13 parts), array_cwise, vectorization_logic, stable_norm, packetmath_1/2, cholesky, realview, block, ref, nomalloc: 125/125 CTest cases pass under GCC 13.3 / AVX2. redux_* additionally passes at r15 across seeds 1, 7, 42, 1234 and 99999, and the whole set passes under GCC / SSE2 and Clang 18.1.3 / AVX2.

test/redux.cpp's redux_vec_boundary covered 4 * PS - 1, 4 * PS, 4 * PS + 1, 8 * PS and 8 * PS + 1, which pin the loop bound but never the residual branches: those need 5 * PS, 6 * PS and 7 * PS. The third commit adds them, and runs every size again through a Map at each offset in 1..PS - 1, since Vec::Random() is always default-aligned and nothing else reaches alignedStart > 0 together with the four-way loop.

Four mutations of the new code -- dropping either of the two- and three-packet residual folds, dropping the third packet of the straight-line path, and starting the four-way loop one packet early -- each fail redux_12 now. Before that commit, the first of them passed five runs at independent seeds.

The sliced path gets redux_block_boundary: inner sizes from one packet to 8 * PS + 1 against one, two, three, four, five, eight and nine panels, for float, double and complex<float>, checking sum, squaredNorm, minCoeff and maxCoeff against scalar references. Four mutations of its three paths fail it as well -- dropping the wide path's third residual, skipping the narrow path's leftover panels, missing the first four panels' packets beyond the first, and starting panel zero one packet early. block, corners and visitor pass alongside the sets above.

Changed-line clang-format 17 is clean.

Note for reviewers

This regroups floating-point reductions from two partial results into four, so sum(), mean() and the norms move in their last bits. Reductions were already reassociated relative to a serial loop, so this is a change of grouping rather than a new guarantee being broken, and the regrouping tightens the bound rather than loosening it: the leading term of the summation error goes from \tfrac{n}{2}\,u\sum|x_i| with two accumulators to \tfrac{n}{4}\,u\sum|x_i| with four. Still worth a deliberate nod before merging. Block reductions go from one partial result to four, where the same bound improves from n\,u\sum|x_i| to \tfrac{n}{4}\,u\sum|x_i|.

I found this while reviewing !3108 (merged), whose GCC double PropagateFast regressions turn out to originate here rather than in that MR's gather work. Details are in a comment there. !3108 (merged) does not depend on this MR, but it benefits from it.

🤖 Generated with Claude Code