Optimization: AVX2 vectorization of the ped-ped loop (float, 8-wide)

Author: Linus
Predecessor: variant_float_inlining_SoA.cpp
File(s): variants/variant_float_inlining_SoA_AVX2.cpp
Status: TBD

What this changes

Float port of the double variant_inlining_7_SoA_AVX2. Inner ped-ped j loop runs on __m256 (8 floats) instead of the double version's __m256d (4 doubles). 8 ped-ped interactions per j step. SoA layout from the predecessor makes the j-side loads contiguous.

Same masking trick as the double variant: build an active mask from j_idx != i AND j_idx < n_p, blend a safe 1.0 into the denominator where active is zero so the i==j and padding lanes don't NaN out, then vandps the per-lane f_x/f_y against active before accumulating into acc_x/acc_y.

n_p is padded to a multiple of 8 (n_p_pad = (n_p + 7) & ~7) so the j loop has no scalar tail. Padding positions are set to 1e9f and sb_e, sb_sq to 0 so any unmasked contribution from those lanes would still be ~zero through expf(very negative) ≈ 0.

EPS clamps from the float scalar predecessor are kept and lifted to SIMD form (_mm256_max_ps(r_ab_n, v_eps) etc.). Same correctness reason as the scalar version: float precision around close encounters can collapse norms to zero.

The expf calls are still scalar inside the j loop: 8 of them per j step, manually unrolled into a scratch tmp_arg[8]/tmp_val[8] array. gcc refused to unroll the loop on its own (it serialised all 8 calls on xmm0), the manual unroll at least lets them issue independently. Kill these next with the polyexp variant.

Code (relevant snippet)

for (int j = 0; j < n_p_pad; j += 8) {
  __m256 j_idx = _mm256_add_ps(_mm256_set1_ps((float)j), lane_off);
  __m256 active = _mm256_and_ps(
      _mm256_cmp_ps(j_idx, v_i,  _CMP_NEQ_OQ),
      _mm256_cmp_ps(j_idx, v_np, _CMP_LT_OQ));

  __m256 pjx    = _mm256_loadu_ps(&pos_x[j]);
  __m256 pjy    = _mm256_loadu_ps(&pos_y[j]);
  // ... sb_e_x, sb_e_y, sb_sq

  // body: r_ab, r_ab2, two_b, scalar expf x 8, smc, sm_r_ab, sm_r_ab2,
  // f_x, f_y, w_val branch, AND with active, fmadd into acc_x/acc_y.
}

See variants/variant_float_inlining_SoA_AVX2.cpp for the full inner loop.

Benchmark

Conditions: min of N runs (TBD), N_ITER=1000, -O3 -march=native -flto. Diff is vs variant_float_inlining_SoA.

x86

Machine: TBD

Frequency: TBD

N	W(n) [Flops]	T(n) [Cycles]	Performance [Flops/Cycle]	Diff vs Predecessor
50
100
200
500

ARM (Apple)

Not applicable, this variant uses AVX2 intrinsics.

Analysis

Same shape of win as the double AVX2 variant, just at 8-wide. Scalar expf is still the dominant cost, polyexp kills it next.

One catch: this variant only hits its real performance if FTZ is enabled in main-float.cpp (that gets added on the polyexp MR). Without it the active-lane chain produces subnormal intermediates (the smc = numer/denom step is the one that lands there for typical mid-range pair distances) and every following FMA stalls 100+ cycles in microcode. With FTZ on it runs as expected, without it benchmarks look 5-10x slow.

Verify passes against baseline_float_n60_i500.csv at 1cm tolerance. Local max drift ~3mm.

Takeaway

Polyexp is the next step. Note: benchmarks are only meaningful with FTZ/DAZ on in main-float.cpp.

3.7 KiB Raw Blame History