-O2                       loop vectorizer + SLP on
-fno-vectorize            disable loop vectorizer
-fno-slp-vectorize        disable straight-line (superword) vectorizer
-march=x86-64-v2/v3/v4    portable ISA levels: SSE4.2 / AVX2+FMA+BMI / AVX-512
-mprefer-vector-width=N   override LLVM's default cap
-mtune=…                  cost model only, no new instructions

for (size_t i = 0; i < n; i++) a[i] = b[i] * k;

-O2                                        width 4,  interleave 2   (SSE2 baseline)
-O2 -march=native                          width 8,  interleave 4   (256-bit, not 512!)
-O2 -march=native -mprefer-vector-width=512  width 16, interleave 4

clang -O2 -march=x86-64-v3 \
      -Rpass=loop-vectorize \         # what vectorized
      -Rpass-missed=loop-vectorize \  # what didn't
      -Rpass-analysis=loop-vectorize  # why not
clang -O2 -fsave-optimization-record -foptimization-record-file=rec.yaml -c f.c

cmpq $32, %rdx      ; n >= 32 ?
jb   .LBB0_6        ; -> scalar path
cmpq $128, %rcx     ; |a-b| far enough apart ?
jb   .LBB0_6

text size, same loop:   no restrict 315 bytes    with restrict 166 bytes

remark: loop not vectorized: cannot prove it is safe to reorder floating-point
        operations; allow reordering by specifying '#pragma clang loop vectorize(enable)'
        before the loop or by providing the compiler option '-ffast-math'

/* (a) per loop — also implies "reorder FP here" */
#pragma clang loop vectorize(enable)
for (size_t i = 0; i < n; i++) s += a[i];

/* (b) per block — MUST start a compound statement, not merely precede a loop */
{
#pragma clang fp reassociate(on)
    for (size_t i = 0; i < n; i++) s += a[i];
}

/* (c) OpenMP semantics, no runtime needed:  clang -fopenmp-simd */
#pragma omp simd reduction(+:s)
for (size_t i = 0; i < n; i++) s += a[i];

/* (d) -ffast-math — whole TU, and it links crtfastmath.o which sets FTZ/DAZ process-wide */

a*b + c   default / =fast :  vfmadd213ps
          -ffp-contract=off: vmulps ; vaddps

plain:                   remark: library call cannot be vectorized. Try compiling
                                 with -fno-math-errno, -ffast-math, or similar flags
-fno-math-errno:         remark: vectorized loop (width 8) — but the asm holds 9 scalar
                                 `callq sinf`: vectorized *around* scalarized calls
-fveclib=libmvec
     + -fno-math-errno:  callq _ZGVdN8v_sinf   ← a real 8-wide vector sine

a = __builtin_assume_aligned(a, 64);
b = __builtin_assume_aligned(b, 64);
n &= ~(size_t)15;                      /* trip count is a multiple of 16 */
for (size_t i = 0; i < n; i++) a[i] += b[i];

#pragma clang loop vectorize(enable|disable)
#pragma clang loop vectorize_width(8)            // or (8, fixed) / (4, scalable)
#pragma clang loop interleave_count(4)           // accumulator count
#pragma clang loop vectorize_predicate(enable)   // fold the tail into a mask (SVE, AVX-512)
#pragma clang loop unroll_count(4)
#pragma clang loop distribute(enable)            // split a loop to break a dependence

#pragma clang attribute push(__attribute__((target("avx2,fma"))), apply_to=function)
__m256 f(__m256 a, __m256 b, __m256 c){ return _mm256_fmadd_ps(a,b,c); }
#pragma clang attribute pop

typedef float f32x8 __attribute__((vector_size(32)));    /* GCC-compatible */
typedef float e32x8 __attribute__((ext_vector_type(8))); /* OpenCL-style   */

e32x8 z = x > y ? x : y;   /* OK -> lane-wise max */
f32x8 z = x > y ? x : y;   /* error: used type '…vector of 8 int' where arithmetic
                              or pointer type is required */

__attribute__((target("avx2,fma")))  float dot_avx2  (…);
__attribute__((target("avx512f")))   float dot_avx512(…);

for (; i < n; i++) s += a[i]*b[i];                 /* scalar epilogue: always correct */

__mmask16 m = (__mmask16)((1u << (n - i)) - 1u);   /* AVX-512 mask: no epilogue at all */
acc = _mm512_fmadd_ps(_mm512_maskz_loadu_ps(m, a+i),
                      _mm512_maskz_loadu_ps(m, b+i), acc);

/* over-read into padding you own — fastest; a page-fault bug if you don't own it */

c0 = vfmaq_f32(c0, vld1q_f32(a+i), vld1q_f32(b+i));  /* -> fmla v0.4s, v4.4s, v2.4s */
float s = vaddvq_f32(c0);                             /* -> faddp / addv             */

svfloat32_t acc = svdup_f32(0.0f);
for (size_t i = 0; i < n; i += svcntw()) {
    svbool_t pg = svwhilelt_b32(i, n);
    acc = svmla_f32_x(pg, acc, svld1_f32(pg,a+i), svld1_f32(pg,b+i));
}
return svaddv_f32(svptrue_b32(), acc);
/* -> whilelt / fmla z0.s, p0/m, z1.s, z2.s / faddv s0, p0, z0.s */

$ clang --target=aarch64-linux-gnu -march=armv8-a+sve -O2 -Rpass=loop-vectorize
remark: vectorized loop (vectorization width: vscale x 4, interleaved count: 2)

__attribute__((target_clones("default","avx2","avx512f")))
float sum(const float *restrict a, size_t n) { … }

0000000000002600 i sum
00000000000022b0 T sum.avx2.0
0000000000002380 T sum.avx512f.1
0000000000002230 T sum.default.2
0000000000002600 W sum.resolver
R_X86_64_IRELATIV  2600

if (__builtin_cpu_supports("avx512f")) return dot_avx512;
if (__builtin_cpu_supports("avx2") && __builtin_cpu_supports("fma")) return dot_avx2;
return dot_scalar;

#define KERNEL(NAME, W)                                                        \
typedef float NAME##_v __attribute__((vector_size(4*(W))));                    \
__attribute__((always_inline)) static inline                                   \
float NAME##_impl(const float *restrict a, const float *restrict b, size_t n){ \
    NAME##_v acc = {0}; size_t i = 0;                                          \
    for (; i + (W) <= n; i += (W)) {                                           \
        NAME##_v va, vb; memcpy(&va,a+i,4*(W)); memcpy(&vb,b+i,4*(W));         \
        acc = va*vb + acc;                                                     \
    }                                                                          \
    float s = 0; for (int k=0;k<(W);k++) s += acc[k];                          \
    for (; i < n; i++) s += a[i]*b[i];                                         \
    return s;                                                                  \
}
KERNEL(k4, 4) KERNEL(k8, 8) KERNEL(k16, 16)

__attribute__((target("avx512f"))) static float dot_512(…){ return k16_impl(a,b,n); }
__attribute__((target("avx2,fma"))) static float dot_256(…){ return  k8_impl(a,b,n); }
static                              float dot_128(…){ return  k4_impl(a,b,n); }

c[i] = a[i] * b[i] + c[i];        // 8× vmulss + 8× vaddss
c[i] = a[i].mul_add(b[i], c[i]);  // 8× vfmadd213ss

let ca = a.chunks_exact(W);
let cb = b.chunks_exact(W);
let (ra, rb) = (ca.remainder(), cb.remainder());   // grab tails BEFORE consuming
let mut acc = [0.0f32; W];
for (x, y) in ca.zip(cb) {
    for i in 0..W { acc[i] += x[i] * y[i]; }
}

RUSTFLAGS="-C target-cpu=native" cargo build --release
rustc -O -C target-cpu=x86-64-v3        # portable ISA levels work here too
rustc -O -C llvm-args=-pass-remarks=loop-vectorize    # LLVM remarks, no -Rpass equivalent

#![feature(portable_simd)]
use std::simd::{f32x8, num::SimdFloat, StdFloat};  // on 1.75: std::simd::SimdFloat
let mut acc = [f32x8::splat(0.0); 4];              // four chains again
acc[i % 4] = vx.mul_add(vy, acc[i % 4]);

rustc -O                          1.69 GFLOP/s   ← mul_add lowered to 8× callq fmaf
rustc -O -C target-cpu=x86-64-v3 30.62 GFLOP/s
rustc -O -C target-cpu=native    31.58 GFLOP/s

#[target_feature(enable = "avx2,fma")]
unsafe fn dot_avx2(a: &[f32], b: &[f32]) -> f32 { /* _mm256_fmadd_ps ×4 */ }

pub fn dot(a: &[f32], b: &[f32]) -> f32 {
    if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
        return unsafe { dot_avx2(a, b) };
    }
    dot_chunked32(a, b)
}

scalar      2.98 GFLOP/s   .iter().zip().map().sum() — cannot reassociate
chunked8    4.59 GFLOP/s
chunked32  14.77 GFLOP/s
avx2       31.87 GFLOP/s   ← full speed from a baseline build
dispatch   31.84 GFLOP/s
mul_add     1.58 GFLOP/s   ← serial fold of true fmas: one 4-cycle chain

#[inline(always)]                                    // ← this is what makes it work
fn dot_kernel<const W: usize>(a: &[f32], b: &[f32]) -> f32 { /* §11 body */ }

#[target_feature(enable = "avx512f")]  unsafe fn dot_512(…) -> f32 { dot_kernel::<64>(a, b) }
#[target_feature(enable = "avx2,fma")] unsafe fn dot_256(…) -> f32 { dot_kernel::<32>(a, b) }
fn dot_128(…) -> f32 { dot_kernel::<16>(a, b) }

sse2        6.08 GFLOP/s
avx2       28.56 GFLOP/s
avx512     24.19 GFLOP/s   ← slower than AVX2: 512-bit downclocking, exactly as in C
dispatch   23.63 GFLOP/s

clang -O2 -march=x86-64-v3 dot_bench.c -o dot_bench && ./dot_bench
clang -O2                  dot_bench.c -o dot_base  && ./dot_base

/* build: clang -O2 -march=x86-64-v3 06_dot.c -o 06_dot */
#include <immintrin.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

/* ---- 1. scalar reference (strictly ordered) ------------------------ */
__attribute__((noinline))
float dot_scalar(const float *restrict a, const float *restrict b, size_t n) {
    float s = 0.0f;
    for (size_t i = 0; i < n; i++) s += a[i] * b[i];
    return s;
}

/* ---- 2. autovectorised: the pragma licenses FP reassociation ------- */
__attribute__((noinline))
float dot_autovec(const float *restrict a, const float *restrict b, size_t n) {
    float s = 0.0f;
#pragma clang loop vectorize(enable) interleave_count(4)
    for (size_t i = 0; i < n; i++) s += a[i] * b[i];
    return s;
}

/* ---- 3. portable generic vectors (no target header) ---------------- */
typedef float f32v8 __attribute__((vector_size(32)));
typedef float f32v4 __attribute__((vector_size(16)));
static inline float hsum8(f32v8 v) {
    f32v4 lo = __builtin_shufflevector(v, v, 0,1,2,3);
    f32v4 hi = __builtin_shufflevector(v, v, 4,5,6,7);
    f32v4 q  = lo + hi;
    return q[0] + q[1] + q[2] + q[3];
}
__attribute__((noinline))
float dot_generic(const float *restrict a, const float *restrict b, size_t n) {
    f32v8 acc = {0};
    size_t i = 0;
    for (; i + 8 <= n; i += 8) {
        f32v8 va, vb;
        memcpy(&va, a + i, 32);            /* unaligned-safe, no aliasing UB */
        memcpy(&vb, b + i, 32);
        acc = va*vb + acc;   /* -ffp-contract=on gives FMA where available */
    }
    float s = hsum8(acc);
    for (; i < n; i++) s += a[i] * b[i];   /* scalar tail */
    return s;
}


/* ---- 3b. portable generic vectors, 4 accumulators ------------------ */
__attribute__((noinline))
float dot_generic4(const float *restrict a, const float *restrict b, size_t n) {
    f32v8 c0 = {0}, c1 = {0}, c2 = {0}, c3 = {0};
    size_t i = 0;
    for (; i + 32 <= n; i += 32) {
        f32v8 a0,a1,a2,a3,b0,b1,b2,b3;
        memcpy(&a0,a+i,32);    memcpy(&b0,b+i,32);
        memcpy(&a1,a+i+8,32);  memcpy(&b1,b+i+8,32);
        memcpy(&a2,a+i+16,32); memcpy(&b2,b+i+16,32);
        memcpy(&a3,a+i+24,32); memcpy(&b3,b+i+24,32);
        c0 = a0*b0 + c0;  c1 = a1*b1 + c1;      /* contracted to FMA by default */
        c2 = a2*b2 + c2;  c3 = a3*b3 + c3;
    }
    f32v8 acc = (c0+c1) + (c2+c3);
    for (; i + 8 <= n; i += 8) { f32v8 va,vb;
        memcpy(&va,a+i,32); memcpy(&vb,b+i,32); acc = va*vb + acc; }
    float s = hsum8(acc);
    for (; i < n; i++) s += a[i] * b[i];
    return s;
}

/* ---- 4. AVX2 intrinsics, 4 accumulators to hide FMA latency -------- */
__attribute__((noinline, target("avx2,fma")))
float dot_avx2(const float *restrict a, const float *restrict b, size_t n) {
    __m256 a0 = _mm256_setzero_ps(), a1 = a0, a2 = a0, a3 = a0;
    size_t i = 0;
    for (; i + 32 <= n; i += 32) {
        a0 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i   ), _mm256_loadu_ps(b+i   ), a0);
        a1 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+ 8), _mm256_loadu_ps(b+i+ 8), a1);
        a2 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+16), _mm256_loadu_ps(b+i+16), a2);
        a3 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i+24), _mm256_loadu_ps(b+i+24), a3);
    }
    for (; i + 8 <= n; i += 8)
        a0 = _mm256_fmadd_ps(_mm256_loadu_ps(a+i), _mm256_loadu_ps(b+i), a0);
    __m256 v = _mm256_add_ps(_mm256_add_ps(a0,a1), _mm256_add_ps(a2,a3));
    __m128 q = _mm_add_ps(_mm256_castps256_ps128(v), _mm256_extractf128_ps(v,1));
    q = _mm_add_ps(q, _mm_movehl_ps(q,q));
    q = _mm_add_ss(q, _mm_shuffle_ps(q,q,1));
    float s = _mm_cvtss_f32(q);
    for (; i < n; i++) s += a[i] * b[i];
    return s;
}

/* ---- 5. AVX-512 with a masked tail: no scalar epilogue ------------- */
__attribute__((noinline, target("avx512f")))
float dot_avx512(const float *restrict a, const float *restrict b, size_t n) {
    __m512 acc = _mm512_setzero_ps();
    size_t i = 0;
    for (; i + 16 <= n; i += 16)
        acc = _mm512_fmadd_ps(_mm512_loadu_ps(a+i), _mm512_loadu_ps(b+i), acc);
    if (i < n) {
        __mmask16 m = (__mmask16)((1u << (n - i)) - 1u);
        acc = _mm512_fmadd_ps(_mm512_maskz_loadu_ps(m, a+i),
                              _mm512_maskz_loadu_ps(m, b+i), acc);
    }
    return _mm512_reduce_add_ps(acc);
}

/* ---- 6. runtime dispatch ------------------------------------------- */
typedef float (*dot_fn)(const float *, const float *, size_t);
static dot_fn pick(void) {
    if (__builtin_cpu_supports("avx512f")) return dot_avx512;
    if (__builtin_cpu_supports("avx2") && __builtin_cpu_supports("fma")) return dot_avx2;
    return dot_scalar;
}

/* ---- harness -------------------------------------------------------- */
static double now(void){ struct timespec t; clock_gettime(CLOCK_MONOTONIC,&t);
                         return t.tv_sec + 1e-9*t.tv_nsec; }
#define N 8003
int main(void) {
    float *a = aligned_alloc(64, ((N*sizeof(float))+63)/64*64);
    float *b = aligned_alloc(64, ((N*sizeof(float))+63)/64*64);
    for (size_t i=0;i<N;i++){ a[i]=(float)((long)(i%17)-8)*0.25f; b[i]=(float)((long)(i%13)-6)*0.5f; }

    struct { const char *nm; dot_fn f; } v[] = {
        {"scalar ", dot_scalar}, {"autovec", dot_autovec}, {"generic ", dot_generic}, {"generic4", dot_generic4},
        {"avx2   ", dot_avx2},   {"avx512 ", dot_avx512},  {"dispatch", pick()} };

    for (unsigned k=0;k<sizeof v/sizeof*v;k++) {
        float r = v[k].f(a,b,N);
        double t0=now(); volatile float sink=0;
        for (int it=0; it<20000; it++) sink += v[k].f(a,b,N);
        double dt=now()-t0;
        printf("%s  result=%12.4f   %6.2f GFLOP/s\n", v[k].nm, r,
               2.0*N*20000/dt/1e9);
    }
    free(a); free(b); return 0;
}

rustc -O dot.rs -o dot && ./dot
rustc -O -C target-cpu=x86-64-v3 dot.rs -o dot3 && ./dot3

use std::hint::black_box;
use std::time::Instant;

// 1. naive: LLVM may NOT reorder these adds -> stays scalar
#[inline(never)]
pub fn dot_scalar(a: &[f32], b: &[f32]) -> f32 {
    a.iter().zip(b).map(|(x, y)| x * y).sum()
}

// 2. explicit lanes: the array accumulator gives 8 independent chains
#[inline(never)]
pub fn dot_chunked(a: &[f32], b: &[f32]) -> f32 {
    const W: usize = 8;
    let mut acc = [0.0f32; W];
    let ca = a.chunks_exact(W);
    let cb = b.chunks_exact(W);
    let (ra, rb) = (ca.remainder(), cb.remainder());
    for (x, y) in ca.zip(cb) {
        for i in 0..W {
            acc[i] += x[i] * y[i];
        }
    }
    let mut s: f32 = acc.iter().sum();
    for (x, y) in ra.iter().zip(rb) {
        s += x * y;
    }
    s
}

// 3. 32 lanes = 4 vectors of 8 -> 4 dependency chains
#[inline(never)]
pub fn dot_chunked4(a: &[f32], b: &[f32]) -> f32 {
    const W: usize = 32;
    let mut acc = [0.0f32; W];
    let ca = a.chunks_exact(W);
    let cb = b.chunks_exact(W);
    let (ra, rb) = (ca.remainder(), cb.remainder());
    for (x, y) in ca.zip(cb) {
        for i in 0..W {
            acc[i] += x[i] * y[i];
        }
    }
    let mut s: f32 = acc.iter().sum();
    for (x, y) in ra.iter().zip(rb) {
        s += x * y;
    }
    s
}

// 4. intrinsics
#[cfg(target_arch = "x86_64")]
use std::arch::x86_64::*;

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2,fma")]
#[inline(never)]
pub unsafe fn dot_avx2(a: &[f32], b: &[f32]) -> f32 {
    let n = a.len().min(b.len());
    let (mut c0, mut c1) = (_mm256_setzero_ps(), _mm256_setzero_ps());
    let (mut c2, mut c3) = (_mm256_setzero_ps(), _mm256_setzero_ps());
    let (pa, pb) = (a.as_ptr(), b.as_ptr());
    let mut i = 0;
    while i + 32 <= n {
        c0 = _mm256_fmadd_ps(_mm256_loadu_ps(pa.add(i)),      _mm256_loadu_ps(pb.add(i)),      c0);
        c1 = _mm256_fmadd_ps(_mm256_loadu_ps(pa.add(i + 8)),  _mm256_loadu_ps(pb.add(i + 8)),  c1);
        c2 = _mm256_fmadd_ps(_mm256_loadu_ps(pa.add(i + 16)), _mm256_loadu_ps(pb.add(i + 16)), c2);
        c3 = _mm256_fmadd_ps(_mm256_loadu_ps(pa.add(i + 24)), _mm256_loadu_ps(pb.add(i + 24)), c3);
        i += 32;
    }
    while i + 8 <= n {
        c0 = _mm256_fmadd_ps(_mm256_loadu_ps(pa.add(i)), _mm256_loadu_ps(pb.add(i)), c0);
        i += 8;
    }
    let v = _mm256_add_ps(_mm256_add_ps(c0, c1), _mm256_add_ps(c2, c3));
    let mut q = _mm_add_ps(_mm256_castps256_ps128(v), _mm256_extractf128_ps(v, 1));
    q = _mm_add_ps(q, _mm_movehl_ps(q, q));
    q = _mm_add_ss(q, _mm_shuffle_ps(q, q, 1));
    let mut s = _mm_cvtss_f32(q);
    while i < n {
        s += a[i] * b[i];
        i += 1;
    }
    s
}

// 5. runtime dispatch
pub fn dot(a: &[f32], b: &[f32]) -> f32 {
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
            return unsafe { dot_avx2(a, b) };
        }
    }
    dot_chunked4(a, b)
}

// 6. the mul_add trap: this is a TRUE fma, not "multiply then add"
#[inline(never)]
pub fn dot_muladd(a: &[f32], b: &[f32]) -> f32 {
    a.iter().zip(b).fold(0.0f32, |s, (x, y)| x.mul_add(*y, s))
}

const N: usize = 8003;
fn bench(name: &str, f: impl Fn(&[f32], &[f32]) -> f32, a: &[f32], b: &[f32]) {
    let r = f(a, b);
    let t = Instant::now();
    let mut sink = 0.0f32;
    for _ in 0..20_000 {
        sink += black_box(f(black_box(a), black_box(b)));
    }
    let dt = t.elapsed().as_secs_f64();
    println!("{name:10} result={r:12.4}  {:6.2} GFLOP/s  (sink {})",
             2.0 * N as f64 * 20_000.0 / dt / 1e9, sink.is_finite());
}
fn main() {
    let a: Vec<f32> = (0..N).map(|i| ((i % 17) as i64 - 8) as f32 * 0.25).collect();
    let b: Vec<f32> = (0..N).map(|i| ((i % 13) as i64 - 6) as f32 * 0.5).collect();
    bench("scalar", dot_scalar, &a, &b);
    bench("chunked8", dot_chunked, &a, &b);
    bench("chunked32", dot_chunked4, &a, &b);
    bench("avx2", |a, b| unsafe { dot_avx2(a, b) }, &a, &b);
    bench("dispatch", dot, &a, &b);
    bench("mul_add", dot_muladd, &a, &b);
}

rustc -O clone.rs -o clone && ./clone      # needs Rust >= 1.89 for the AVX-512 arm

// Requires Rust >= 1.89 (AVX-512 target features stabilised). Build: rustc -O clone.rs

use std::hint::black_box;
use std::time::Instant;

/// One generic kernel. `#[inline(always)]` is what makes it get *re-codegen'd*
/// with the caller's target features.
#[inline(always)]
fn dot_kernel<const W: usize>(a: &[f32], b: &[f32]) -> f32 {
    let ca = a.chunks_exact(W);
    let cb = b.chunks_exact(W);
    let (ra, rb) = (ca.remainder(), cb.remainder());
    let mut acc = [0.0f32; W];
    for (x, y) in ca.zip(cb) {
        for i in 0..W { acc[i] += x[i] * y[i]; }
    }
    let mut s = 0.0f32;
    for v in acc { s += v; }
    for (x, y) in ra.iter().zip(rb) { s += x * y; }
    s
}

#[target_feature(enable = "avx512f")]
unsafe fn dot_512(a: &[f32], b: &[f32]) -> f32 { dot_kernel::<64>(a, b) }
#[target_feature(enable = "avx2,fma")]
unsafe fn dot_256(a: &[f32], b: &[f32]) -> f32 { dot_kernel::<32>(a, b) }
fn dot_128(a: &[f32], b: &[f32]) -> f32 { dot_kernel::<16>(a, b) }

pub fn dot(a: &[f32], b: &[f32]) -> f32 {
    #[cfg(target_arch = "x86_64")]
    {
        if is_x86_feature_detected!("avx512f") { return unsafe { dot_512(a, b) }; }
        if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
            return unsafe { dot_256(a, b) };
        }
    }
    dot_128(a, b)
}

const N: usize = 8003;
fn bench(name: &str, f: impl Fn(&[f32], &[f32]) -> f32, a: &[f32], b: &[f32]) {
    let r = f(a, b);
    let t = Instant::now(); let mut s = 0.0f32;
    for _ in 0..20_000 { s += black_box(f(black_box(a), black_box(b))); }
    println!("{name:9} result={r:12.4}  {:6.2} GFLOP/s  {}",
        2.0*N as f64*20_000.0/t.elapsed().as_secs_f64()/1e9, s.is_finite());
}
fn main() {
    let a: Vec<f32> = (0..N).map(|i| ((i % 17) as i64 - 8) as f32 * 0.25).collect();
    let b: Vec<f32> = (0..N).map(|i| ((i % 13) as i64 - 6) as f32 * 0.5).collect();
    bench("sse2",     dot_128, &a, &b);
    bench("avx2",     |a,b| unsafe { dot_256(a,b) }, &a, &b);
    bench("avx512",   |a,b| unsafe { dot_512(a,b) }, &a, &b);
    bench("dispatch", dot, &a, &b);
}

rustc -O psimd.rs -o psimd && ./psimd                            #  1.69 GFLOP/s
rustc -O -C target-cpu=x86-64-v3 psimd.rs -o psimd3 && ./psimd3  # 30.62 GFLOP/s

#![feature(portable_simd)]
use std::simd::{f32x8, SimdFloat, StdFloat};
use std::hint::black_box;
use std::time::Instant;

#[inline(never)]
pub fn dot_psimd(a: &[f32], b: &[f32]) -> f32 {
    const W: usize = 8;
    let ca = a.chunks_exact(W);
    let cb = b.chunks_exact(W);
    let (ra, rb) = (ca.remainder(), cb.remainder());
    let mut acc = [f32x8::splat(0.0); 4];
    for (i, (x, y)) in ca.zip(cb).enumerate() {
        let (vx, vy) = (f32x8::from_slice(x), f32x8::from_slice(y));
        acc[i % 4] = vx.mul_add(vy, acc[i % 4]);
    }
    let v = (acc[0] + acc[1]) + (acc[2] + acc[3]);
    let mut s = v.reduce_sum();
    for (x, y) in ra.iter().zip(rb) { s += x * y; }
    s
}
const N: usize = 8003;
fn main() {
    let a: Vec<f32> = (0..N).map(|i| ((i % 17) as i64 - 8) as f32 * 0.25).collect();
    let b: Vec<f32> = (0..N).map(|i| ((i % 13) as i64 - 6) as f32 * 0.5).collect();
    let r = dot_psimd(&a, &b);
    let t = Instant::now(); let mut s = 0.0f32;
    for _ in 0..20_000 { s += black_box(dot_psimd(black_box(&a), black_box(&b))); }
    println!("psimd    result={r:12.4}  {:6.2} GFLOP/s  {}",
        2.0*N as f64*20_000.0/t.elapsed().as_secs_f64()/1e9, s.is_finite());
}

rustc -O -C target-cpu=x86-64-v3 --emit asm contract.rs
  ->  8 vaddss   8 vfmadd213ss   8 vmulss

#[inline(never)]
pub fn madd(a: &[f32], b: &[f32], c: &mut [f32]) {
    for i in 0..c.len() { c[i] = a[i] * b[i] + c[i]; }
}
#[inline(never)]
pub fn madd_fma(a: &[f32], b: &[f32], c: &mut [f32]) {
    for i in 0..c.len() { c[i] = a[i].mul_add(b[i], c[i]); }
}
fn main(){ let a=[1.0f32;8]; let b=[2.0f32;8]; let mut c=[3.0f32;8];
  madd(&a,&b,&mut c); madd_fma(&a,&b,&mut c); println!("{}", c[0]); }

#include <stddef.h>
float sum_strict(const float *a, size_t n) {
    float s = 0.0f;
    for (size_t i = 0; i < n; i++) s += a[i];
    return s;
}
float sum_pragma(const float *a, size_t n) {
    float s = 0.0f;
    {
#pragma clang fp reassociate(on)
        for (size_t i = 0; i < n; i++) s += a[i];
    }
    return s;
}
float sum_omp(const float *a, size_t n) {
    float s = 0.0f;
#pragma omp simd reduction(+:s)
    for (size_t i = 0; i < n; i++) s += a[i];
    return s;
}

#include <stdio.h>
#include <stdint.h>
typedef float  f32x8 __attribute__((vector_size(32)));
typedef int32_t i32x8 __attribute__((vector_size(32)));
typedef float  e32x8 __attribute__((ext_vector_type(8)));
int main(void){
  f32x8 a={1,2,3,4,5,6,7,8}, b={8,7,6,5,4,3,2,1};
  i32x8 ia={1,2,3,4,5,6,7,8};
  printf("reduce_add int   = %d\n", __builtin_reduce_add(ia));
  printf("reduce_max float = %.1f\n", __builtin_reduce_max(a));
  printf("elementwise_max0 = %.1f\n", __builtin_elementwise_max(a,b)[0]);
  e32x8 x={1,2,3,4,5,6,7,8}, y={8,7,6,5,4,3,2,1};
  e32x8 z = x > y ? x : y;                 /* ext_vector_type allows ?: */
  printf("ext ternary lane0=%.1f lane7=%.1f\n", z[0], z[7]);
  printf("ext .hi/.lo? swizzle x.s0=%.1f\n", x.s0);
  return 0;
}

#include <stdio.h>
typedef float f32x8 __attribute__((vector_size(32)));
typedef float f32x4 __attribute__((vector_size(16)));
typedef float f32x2 __attribute__((vector_size(8)));

static inline float hsum_f32x8(f32x8 v) {
    f32x4 lo = __builtin_shufflevector(v, v, 0,1,2,3);
    f32x4 hi = __builtin_shufflevector(v, v, 4,5,6,7);
    f32x4 q  = lo + hi;
    f32x2 a  = __builtin_shufflevector(q, q, 0,1);
    f32x2 b  = __builtin_shufflevector(q, q, 2,3);
    f32x2 d  = a + b;
    return d[0] + d[1];
}
int main(void){ f32x8 v={1,2,3,4,5,6,7,8}; printf("%.1f\n", hsum_f32x8(v)); }