SIMD on the Laptop: NEON Intrinsics, Auto-Vectorisation and the Loop That Would Not Vectorise

Every core in a modern laptop has vector registers that hold four single-precision floats and instructions that operate on all four lanes at once. On AArch64 the extension is NEON (Advanced SIMD): 32 registers of 128 bits. There are three ways to get vector code: trust the compiler, steer it with pragmas and qualifiers, or write intrinsics by hand. This post takes two kernels through all three, shows the assembly clang produced, and measures each version with its data in L1, in L2 and in DRAM on one core of an Apple-silicon laptop (machdep.cpu.brand_string: Apple M1, hw.ncpu: 8). The compiler is Apple clang 17.0.0 (clang-1700.0.13.5) with -O3 -mcpu=native.

Two kernels, three implementations

The kernels are the two simplest in linear algebra. The dot product $s=\sum _i{x}_i{y}_i$ is a reduction; saxpy, $y_i\leftarrow a{x}_i+y_i$, is a streaming update in which no iteration depends on another. That difference turns out to decide almost everything below.

• Scalar
• Auto-vectorised
• NEON intrinsics

• float dot_scalar(const float *x, const float *y, size_t n) {
      float s = 0.f;
  #pragma clang loop vectorize(disable) interleave(disable)
      for (size_t i = 0; i < n; i++) s += x[i] * y[i];
      return s;
  }
  

  The pragma stops clang from vectorising the baseline. Every iteration's fmadd depends on the previous one through s.

• /* -O3 alone: a floating-point reduction may not be reassociated */
  float dot_auto(const float *x, const float *y, size_t n) {
      float s = 0.f;
      for (size_t i = 0; i < n; i++) s += x[i] * y[i];
      return s;
  }
  
  /* the same loop, with reassociation allowed inside this block only */
  float dot_reassoc(const float *x, const float *y, size_t n) {
      float s = 0.f;
      {
  #pragma clang fp reassociate(on)
          for (size_t i = 0; i < n; i++) s += x[i] * y[i];
      }
      return s;
  }
  

  A vectorised reduction keeps one partial sum per lane and combines them at the end, which is a different order of additions from the serial loop and therefore, in floating point, potentially a different answer. The second version gives clang permission to do that for this one loop.

• #include <arm_neon.h>
  
  float dot_neon(const float *x, const float *y, size_t n) {
      float32x4_t a0 = vdupq_n_f32(0.f), a1 = a0, a2 = a0, a3 = a0;
      size_t i = 0;
      for (; i + 16 <= n; i += 16) {
          float32x4x4_t xv = vld1q_f32_x4(x + i), yv = vld1q_f32_x4(y + i);   /* 64 bytes each */
          a0 = vfmaq_f32(a0, xv.val[0], yv.val[0]);
          a1 = vfmaq_f32(a1, xv.val[1], yv.val[1]);
          a2 = vfmaq_f32(a2, xv.val[2], yv.val[2]);
          a3 = vfmaq_f32(a3, xv.val[3], yv.val[3]);
      }
      float s = vaddvq_f32(vaddq_f32(vaddq_f32(a0, a1), vaddq_f32(a2, a3)));
  #pragma clang loop vectorize(disable) interleave(disable)
      for (; i < n; i++) s += x[i] * y[i];                  /* tail: n mod 16 elements */
      return s;
  }
  

  Four accumulators, so four multiply-adds can be in flight instead of waiting on each other; vaddvq_f32 adds across the lanes at the end, and a scalar loop takes the last n mod 16 elements.

Saxpy gets the same treatment: a scalar baseline with vectorisation disabled, the plain loop at -O3 with and without restrict, and intrinsics versions that process one vector or four vectors per iteration. The four-vector version is below; I also measured it with twelve separate vld1q_f32/vst1q_f32 calls (eight loads, four stores) in place of the _x4 forms, for a reason that shows up in the measurements.

void saxpy_neon4(float a, const float *restrict x, float *restrict y, size_t n) {
    size_t i = 0;
    for (; i + 16 <= n; i += 16) {
        float32x4x4_t xv = vld1q_f32_x4(x + i), yv = vld1q_f32_x4(y + i);
        yv.val[0] = vfmaq_n_f32(yv.val[0], xv.val[0], a);
        yv.val[1] = vfmaq_n_f32(yv.val[1], xv.val[1], a);
        yv.val[2] = vfmaq_n_f32(yv.val[2], xv.val[2], a);
        yv.val[3] = vfmaq_n_f32(yv.val[3], xv.val[3], a);
        vst1q_f32_x4(y + i, yv);
    }
#pragma clang loop vectorize(disable) interleave(disable)
    for (; i < n; i++) y[i] = a * x[i] + y[i];            /* tail */
}

What clang emitted

The compiler reports its decisions when asked: -Rpass=loop-vectorize -Rpass-missed=loop-vectorize -Rpass-analysis=loop-vectorize. For a file containing the four dot-product functions above followed by saxpy with and without restrict, it printed (remarks about the intrinsics function's own loops omitted):

simd.c:7:5: remark: loop not vectorized: vectorization and interleaving are explicitly disabled, or the loop has already been vectorized [-Rpass-analysis=loop-vectorize]
simd.c:13:5: remark: vectorized loop (vectorization width: 4, interleaved count: 4) [-Rpass=loop-vectorize]
simd.c:21:9: remark: vectorized loop (vectorization width: 4, interleaved count: 4) [-Rpass=loop-vectorize]
simd.c:43:5: remark: vectorized loop (vectorization width: 4, interleaved count: 4) [-Rpass=loop-vectorize]
simd.c:47:5: remark: vectorized loop (vectorization width: 4, interleaved count: 4) [-Rpass=loop-vectorize]

Line 13 is dot_auto, and the remark says it was vectorised. The assembly (clang -O3 -S) shows what that means for a reduction that must keep its order:

LBB1_7:
	ldp	q1, q2, [x9, #-32]
	ldp	q3, q4, [x9], #64
	ldp	q5, q6, [x10, #-32]
	ldp	q7, q16, [x10], #64
	fmul.4s	v1, v1, v5          ; four products at once ...
	mov	s5, v1[3]
	mov	s17, v1[2]
	mov	s18, v1[1]
	...
	fadd	s0, s0, s1          ; ... then sixteen additions, one after another
	fadd	s0, s0, s18
	fadd	s0, s0, s17
	fadd	s0, s0, s5
	...
	subs	x11, x11, #16
	b.ne	LBB1_7

The multiplications are vectorised; the products are then pulled out lane by lane and added one at a time into s0, because the program asked for the exact left-to-right sum. With reassociation allowed, the loop becomes four fmla.4s into four accumulators, which is what the intrinsics version does by hand. Its inner loop shows the instruction the intrinsic is named after:

LBB3_3:
	ld1.4s	{ v4, v5, v6, v7 }, [x9], #64      ; vld1q_f32_x4(x + i)
	ld1.4s	{ v16, v17, v18, v19 }, [x8], #64  ; vld1q_f32_x4(y + i)
	fmla.4s	v0, v16, v4
	fmla.4s	v1, v17, v5
	fmla.4s	v2, v18, v6
	fmla.4s	v3, v19, v7
	add	x11, x10, #32
	add	x10, x10, #16
	cmp	x11, x2
	b.ls	LBB3_3

With the single-register vld1q_f32, clang emits ldr q or a paired ldp q instead of ld1: for one register they load the same bytes. Saxpy with restrict is the textbook case, one source loop becoming three:

LBB4_6:                         ; 16 floats per iteration
	ldp	q1, q2, [x9, #-32]
	ldp	q3, q4, [x9], #64
	ldp	q5, q6, [x10, #-32]
	ldp	q7, q16, [x10]
	fmla.4s	v5, v1, v0[0]
	fmla.4s	v6, v2, v0[0]
	fmla.4s	v7, v3, v0[0]
	fmla.4s	v16, v4, v0[0]
	stp	q5, q6, [x10, #-32]
	stp	q7, q16, [x10], #64
	subs	x11, x11, #16
	b.ne	LBB4_6
LBB4_10:                        ; 4 floats per iteration: vector remainder
	ldr	q1, [x11], #16
	ldr	q2, [x10]
	fmla.4s	v2, v1, v0[0]
	str	q2, [x10], #16
	adds	x9, x9, #4
	b.ne	LBB4_10
LBB4_13:                        ; 1 float per iteration: scalar tail
	ldr	s1, [x10], #4
	ldr	s2, [x8]
	fmadd	s1, s0, s1, s2
	str	s1, [x8], #4
	subs	x9, x9, #1
	b.ne	LBB4_13

Without restrict, clang still vectorises (the remark on line 47 is identical) but first compares the two address ranges and branches to a scalar loop if they overlap.

Measurements

Each timing processes $2^{28}$ elements (many repetitions of a small array or a few passes over a large one), with the arrays sized to sit in L1 (16 KiB each), L2 (1 MiB each) or DRAM (64 MiB each). All kernels are timed round-robin within each round, eleven rounds per pass, three passes; the table gives the median of the three per-pass medians.

Reading down the columns:

• The scalar dot product is latency-bound. 1.54 GFLOP/s matches the 1.52 GFLOP/s of a single dependent scalar FMA chain in the roofline post, and it is the same from L1, L2 and DRAM because the loads were never the bottleneck.
• -O3 alone bought 36 %. That is what remains when the serial chain consists of additions instead of fused multiply-adds.
• Permission to reassociate bought 17 times. Four accumulators of four lanes each is sixteen independent sums. The hand-written intrinsics give the same number: once clang may reassociate, it produces the loop the intrinsics spell out by hand.
• Four accumulators are still a limit in L1. A version with eight accumulators ran at 31.6 GFLOP/s in L1 against 26.7 for four in the same run, and at 20.0 against 20.1 in L2, where bandwidth rather than latency sets the pace.
• Saxpy vectorises on its own, 4.2 times faster than scalar in L1, because nothing in the loop carries from one iteration to the next. Dropping restrict cost nothing measurable: the overlap check runs once per call.
• One vector per iteration is not enough, and the instruction form matters. The single-vector intrinsics loop spends nine instructions on every four floats and reaches 20.8 GFLOP/s; the compiler spends twelve on sixteen. The four-vector version with four-register ld1/st1 reached 24.0, and the identical loop written with single-register loads, which clang pairs into ldp/stp, reached 31.1, the same as the compiler. For the dot product, which only loads, the two forms tied. I have no hardware counters on this machine to say which part of the four-register form costs more; the measurement is that it did.
• From DRAM the code stops mattering. Every vectorised saxpy lands at 6.6 to 6.8 GFLOP/s, about 40 GB/s at 12 bytes per element, and every vectorised dot product near 10 GFLOP/s, also about 40 GB/s. That is about what the triad measured from DRAM in the roofline post (medians of 37 to 40 GB/s). The scalar saxpy, at 4.7 GFLOP/s or 28 GB/s, evidently cannot keep enough memory traffic in flight to reach that bandwidth, so vectorisation still buys 1.4 times there.

WARNING

-ffast-math would also have vectorised the dot product, but it is a bundle: reassociation, the assumption that no NaN or infinity ever appears, flushing denormals and more, applied to the whole translation unit. The reassociated sum is also a different number. On the 4 096-element arrays above the serial loop returned 2170.902 and the reassociated loop 2170.850, against 2170.850 accumulated in double precision; here the reordered sum happened to be the more accurate one, but the point is that the answer changed. #pragma clang fp reassociate(on) scopes the permission to the loop that needs it.

Alignment and tails

NEON loads and stores do not require alignment on normal memory, and misaligned accesses are legal. To see what they cost I offset both arrays by 4 bytes (aligned only to the element) and by 16 bytes (aligned to the register, not to the 128-byte cache line) from a 16 KiB-aligned start:

Only one case stood out: saxpy in L1 lost 27 % when its arrays were offset by 4 bytes. The dot product, which never stores, did not care, and from DRAM nothing rose above the run-to-run noise. Aligned allocation (aligned_alloc(64, ...) or larger) is cheap insurance, and when a function cannot see how its arrays were allocated, __builtin_assume_aligned passes the knowledge along.

The tail is the part people forget. NEON has no masked loads or stores (SVE and AVX-512 do), so a loop over n elements becomes a vector loop over most of them followed by scalar code for the rest, which is exactly the three-loop structure clang generated for saxpy. For long arrays the tail is noise; for a dot product of seven elements it is the whole computation, and the only cure is to arrange data so that the hot loops see long runs.

When vectorisation fails

The remarks also say why a loop stayed scalar. Three loops that do, with clang's own words:

void prefix_sum(float *restrict a, const float *restrict b, size_t n) {
    for (size_t i = 1; i < n; i++)
        a[i] = a[i - 1] + b[i];
}
void prefix_sum_reg(float *restrict a, const float *restrict b, size_t n) {
    float s = 0.f;
    for (size_t i = 0; i < n; i++) { s += b[i]; a[i] = s; }
}
void scatter_add(float *a, const float *b, const int *idx, size_t n) {
    for (size_t i = 0; i < n; i++)
        a[idx[i]] += b[i];
}

fail.c:5:14: remark: loop not vectorized: unsafe dependent memory operations in loop. Use #pragma clang loop distribute(enable) to allow loop distribution to attempt to isolate the offending operations into a separate loop
Backward loop carried data dependence. Memory location is the same as accessed at fail.c:5:16 [-Rpass-analysis=loop-vectorize]
fail.c:4:5: remark: loop not vectorized [-Rpass-missed=loop-vectorize]
fail.c:10:5: remark: loop not vectorized: value that could not be identified as reduction is used outside the loop [-Rpass-analysis=loop-vectorize]
fail.c:10:5: remark: loop not vectorized [-Rpass-missed=loop-vectorize]
fail.c:15:9: remark: loop not vectorized: cannot identify array bounds [-Rpass-analysis=loop-vectorize]
fail.c:14:5: remark: loop not vectorized [-Rpass-missed=loop-vectorize]

-Rpass-missed only says that a loop was not vectorised; -Rpass-analysis says why, so use both. The two prefix sums are the same loop-carried dependency, once through memory and once through a register: each element needs the previous result, so there is nothing to do four at a time without changing the algorithm (a parallel scan does it with extra work). The scatter is the aliasing case. Two indices may be equal, so two lanes could update the same element, and no restrict qualifier can rule that out because the collision happens inside one array; the remark is that clang cannot even bound which addresses the loop touches. The reasons collected in this post:

Further reading

• Arm, Intrinsics reference: every v*q_f32 function and the instruction it maps to.
• Auto-Vectorization in LLVM, on the loop and SLP vectorisers, reductions and runtime alias checks seen above.
• Clang language extensions for #pragma clang loop and #pragma clang fp, and the Clang user's manual for the -Rpass family of optimisation remarks.
• Arm optimized-routines: production NEON code for string and math routines, a good source of idioms for tails and alignment.
• Sergey Slotin, SIMD Parallelism in Algorithms for Modern Hardware, the same ground on x86.