<p>If you need precision <em>and</em> parallelism, use Kahan summation or another error-compensation technique to let you reorder your sum (into SIMD vector element strides with multiple accumulators).</p>
<p>As <a href="https://arxiv.org/pdf/1401.0248.pdf" rel="nofollow noreferrer">Twofold fast summation - Evgeny Latkin</a> points out, if you bottleneck on memory bandwidth, an error-compensated sum isn't much slower than a max-performance sum, since the CPU has lots of computation throughput that goes unused in a simply-parallelized sum that bottlenecks on memory bandwidth</p>
<p>See also (google results for <code>kahan summation avx</code>)</p>
<ul>
<li><p><a href="https://github.com/rreusser/summation-algorithms" rel="nofollow noreferrer">https://github.com/rreusser/summation-algorithms</a></p>
</li>
<li><p><a href="https://scicomp.stackexchange.com/questions/10869/which-algorithm-is-more-accurate-for-computing-the-sum-of-a-sorted-array-of-numb">https://scicomp.stackexchange.com/questions/10869/which-algorithm-is-more-accurate-for-computing-the-sum-of-a-sorted-array-of-numb</a></p>
</li>
</ul>
<hr />
<p>Re: your idea: Summing groups of 4 numbers in-order will let you hide the FP-add latency, and bottleneck on scalar add throughput.</p>
<p>Doing horizontal sums within vectors takes a lot of shuffling, so it's a potential bottleneck.  You might consider loading <code>a0 a1 a2 a3</code>, then shuffling to get <code>a0+a1 x a2+a3 x</code>, then <code>(a0+a1) + (a2+a3)</code>.  You have a Ryzen, right?  The last step is just a <code>vextractf128</code> down to 128b.  That's still 3 total ADD uops, and 3 shuffle uops, but with fewer instructions than scalar or 128b vectors.</p>
<hr />
<h3>Your idea is very similar to Pairwise Summation</h3>
<p>You're always going to get <em>some</em> rounding error, but adding numbers of similar magnitude minimizes it.</p>
<p><strong>See also</strong> <a href="https://stackoverflow.com/questions/55477701/simd-matmul-program-gives-different-numerical-results">Simd matmul program gives different numerical results</a> for some comments about Pairwise Summation and simple efficient SIMD.</p>
<p>The difference between adding 4 contiguous numbers vs. vertically adding 4 SIMD vectors is probably negligible.  SIMD vectors give you small strides (of SIMD vector width) in the array.  Unless the array grows extremely quickly, they're still going to have basically similar magnitudes.</p>
<p><strong>You don't need to horizontal sum until the very end to still get most of the benefit</strong>.  You can maintain 1 or 2 SIMD vector accumulators while you use more SIMD registers to sum short runs (of maybe 4 or 8 SIMD vectors) before adding into the main accumulators.</p>
<p>In fact having your total split more ways (across the SIMD vector elements) means it doesn't grow as large.  So it helps with exactly the problem you're trying to avoid, and <strong>horizontal summing down to a single scalar accumulator actually makes things worse, especially for a strictly sorted array.</strong></p>
<p>With out-of-order execution, you don't need very many tmp accumulators to make this work and hide the FP-add latency of accumulating into the main accumulators.  You can do a couple groups of accumulating into a fresh <code>tmp = _mm_load_ps()</code> vector and adding that to the total, and OoO exec will overlap those executions.  So you don't need a huge unroll factor for your main loop.</p>
<p>But it shouldn't be too small, you don't want to bottleneck on the add into the main accumulator.  You want to bottleneck on FP-add throughput.  (Or if you care about Broadwell/Haswell and you don't totally bottleneck on memory bandwidth, mix in some FMA with a <code>1.0</code> multiplier to take advantage of that throughput.)</p>
<p>e.g. Skylake SIMD FP add has 4 cycle latency, 0.5 cycle throughput, so you need to be doing at least 7 adds that are part of a short dep chain for every add into a single accumulator.  Preferably more, and/or preferably with 2 long-term accumulators to better absorb bubbles in scheduling from resource conflicts.</p>
<p>See <a href="https://stackoverflow.com/questions/66260651/mm256-fmadd-ps-is-slower-than-mm256-mul-ps-mm256-add-ps">_mm256_fmadd_ps is slower than _mm256_mul_ps + _mm256_add_ps?</a> for more about multiple accumulators.</p>
Efficient stable sum of a sorted array in AVX2
<p>Consider a sorted (ascending) array of <code>double</code> numbers. For numerical stability the array should be summed up as if iterating it from the beginning till the end, accumulating the sum in some variable.</p>

<p>How to vectorize this efficiently with AVX2?</p>

<p>I've looked into this method <a href="https://stackoverflow.com/questions/9775538/fastest-way-to-do-horizontal-vector-sum-with-avx-instructions">Fastest way to do horizontal vector sum with AVX instructions</a> , but it seems quite tricky to scale it to an array (some divide&amp;conquer approach may be needed), while keeping the floating-point precision by ensuring that small numbers are summed up before adding them to a larger number.</p>

<p><strong>Clarification 1</strong>:  I think it should be ok to e.g. sum the first 4 items, then add them to the sum of the next 4 items, etc. I'm willing to trade some stability for performance. But I would prefer a method that doesn't ruin the stability completely.</p>

<p><strong>Clarification 2</strong>: memory shouldn't be a bottleneck because the array is in L3 cache (but not in L1/L2 cache, because pieces of the array were populated from different threads). I wouldn't like to resort to Kahan summation because I think it's really the number of operations that matters, and Kahan summation would increase it about 4 times.</p>