Make patch 0005 near upstream speed with in-line F16 checks (#174)

docs/linux_arm64_cpu_variants.md

@@ -35,9 +35,9 @@ enables. The compiler also uses these where they apply:
   Neither is a scored feature; ARMv8.6 guarantees both, and the I8MM that
   scores this variant is an ARMv8.6 feature.
 - `linux_armv8.6_1` and `android_armv8.6_1` also use FMLAL/FMLAL2 (FHM) in the
-  F16 dot kernel of carried patch 0005: FP16 on ARMv8.4 or later implies FHM,
-  so the compiler defines `__ARM_FEATURE_FP16_FML`. The `armv8.2_2` variants
-  use the conversion path instead.
+  F32 recompute of carried patch 0005's F16 dot kernel: FP16 on ARMv8.4 or
+  later implies FHM, so the compiler defines `__ARM_FEATURE_FP16_FML`. The
+  `armv8.2_2` variants convert to F32 there instead.
 - `libggml-cpu.so` (ARMv8.0) contains LSE only inside libgcc's outline-atomic
   helpers, each of which tests `__aarch64_have_lse_atomics` before reaching the
   LSE instruction. It contains no LRCPC, LRCPC2 or CRC.

docs/native174_f16_dot_accumulation.md

@@ -33,10 +33,14 @@ public-model failure is established.
 Both patches change build-tree copies only.
 
 - 0005: under NEON, FMA, IEEE F16, FP16 vector arithmetic and no SVE,
-  `ggml_vec_dot_f16` multiplies, accumulates and reduces in F32 with four
-  low/high accumulator pairs (the original 32-element unroll), an eight-element
-  tail and a scalar tail. With `__ARM_FEATURE_FP16_FML` it uses FMLAL/FMLAL2;
-  otherwise it converts and uses F32 FMA.
+  `ggml_vec_dot_f16` keeps upstream's code for rows of one block (at most
+  2048 elements, so 64 products per lane): F16 multiply-adds in four
+  `float16x8_t` accumulators, the F16 reduction, then F32 and a double tail.
+  If the reduced F16 sums hold an infinity or NaN, the row is computed again
+  with F32 products and sums (FMLAL/FMLAL2 with `__ARM_FEATURE_FP16_FML`,
+  conversion and F32 FMA otherwise). Longer rows run out of line: each block
+  is summed the same way and added into F32, and a nonfinite F32 sum is
+  recomputed the same way.
 - 0006: the FP16 llamafile tile keeps its F16 multiply-adds, but each lane
   sums at most 64 products (512 K elements) before the sum is added into an
   F32 accumulator (FMLAL/FMLAL2 by one with FP16 FML, conversion and F32 add
@@ -45,6 +49,20 @@ Both patches change build-tree copies only.
   `tinyBLAS<4, float32x4_t, ...>` (F16 loads, F32 accumulation), which
   overwrites it.
 
+### Dot contract (0005)
+
+- Overflow: for finite F16 inputs no result is nonfinite. A row whose F16
+  sums overflowed, or with a NaN or infinite input before its scalar tail,
+  comes from the F32 recompute and is bit-identical to the first version of
+  0005. NaN and infinite inputs still give NaN or infinity; in the scalar
+  tail they propagate as upstream.
+- Precision: a finite result for `n` < 2080 (one block plus the scalar tail)
+  is bit-identical to unpatched upstream. Longer rows add each block's F16
+  sums into F32, so F16 rounding no longer compounds across blocks; within a
+  block it is upstream's, so results are not F32-accurate.
+- Cost: a finiteness check of the reduced F16 vector per row; a recomputed
+  row costs the F16 pass plus the F32 pass.
+
 ### GEMM contract (0006)
 
 - Overflow: for finite F16 inputs no output is nonfinite. A tile that
@@ -65,26 +83,38 @@ maintainer approved the
 
 ## Performance
 
-Apple M4 Max, CPU only, `GGML_NATIVE=OFF`. Kernel costs are against the
-unpatched kernels; model throughput is against the shipped `v0.8.0-3` macOS
-library (Qwen2.5-0.5B-Instruct, 8 threads, F16 KV, median of 3 interleaved
-rounds' ratios):
+Model throughput, CPU only: Qwen2.5-0.5B-Instruct, 8 threads, F16 KV, ratio
+to the same llama.cpp `v0.9.0` build without 0005 and 0006, median over
+interleaved rounds (5 on Apple M4 Max, 3 on Google Cloud `t2a` and `c4a`):
+
+| Host, module | F16 pp512, FA off / on | F16 tg128, tg64 at depth 1024 | Q8_0 pp512, FA off / on | Q8_0 tg |
+| --- | --- | --- | --- | --- |
+| M4 Max, macOS arm64 | 0.92 / 0.94 | 0.99–1.00 | 0.98 / 1.00 | 0.99–1.00 |
+| Neoverse-N1, `linux_armv8.2_2` | 1.47 / 1.47 | 0.98–0.99 | 1.05 / 1.00 | 0.99–1.01 |
+| Neoverse-V2, `linux_armv8.2_2` | 1.56 / 1.54 | 0.98–1.00 | 1.09 / 1.00 | 0.99–1.01 |
+| Neoverse-V2, `linux_armv8.6_1` | 0.98 / 0.99 | 0.99–1.02 | 0.98 / 1.00 | 1.00–1.01 |
+
+0006 alone gives the same F16 pp512 on M4 and `linux_armv8.2_2`; 0005 costs
+about 2% of `linux_armv8.6_1` prefill, where llamafile is compiled out and
+every F16 matmul runs `ggml_vec_dot_f16`. The first version of 0005 ran at
+0.50–0.53 for F16 decode on `linux_armv8.2_2` and 0.75 for F16 prefill on
+`linux_armv8.6_1`.
+
+Kernels on M4 Max, one thread, time against the unpatched kernels:
 
 | Case | Cost |
 | --- | --- |
-| F16 mat-vec 4096x4096, 1 thread (0005): FMLAL / conversion path | +13% / +78% time |
-| F16 GEMM 4096x4096, 32–512 columns, 1 or 8 threads (0006): FMLAL / conversion flush | +8–13% / +12–14% time |
-| Same with every tile overflowing (all redone in F32) | +257% time |
-| F16-weight model pp512, FA off / on | 0.91 / 0.94 tokens/s |
-| F16-weight model tg128 and tg64 at depth 1024 (decode does not reach the tile) | 0.96–1.02 (noise) |
-| Q8_0 model pp512, FA off / on | 0.98 / 1.01 tokens/s |
-| Qwen3-TTS cancel latency (core `tts` pack, CPU), median / p90 / max | 16 / 40 / 49 ms against 15 / 32 / 48 ms |
+| F16 mat-vec, rows of 64 / 128 / 256 / 4096 elements (0005): FMLAL / conversion build | +4% / +12% / +2% / −3%; +2% / +5% / +1% / −3% |
+| Same, 4096x4096 with every row overflowing (all recomputed) | 2.1 / 2.75 times |
+| First 0005, 4096x4096 | +14% / +78% |
+| F16 GEMM 4096x4096, 32–512 columns, 1 or 8 threads (0006): FMLAL / conversion flush | +8–13% / +12–14% |
+| Same with every tile overflowing (all redone in F32) | +257% |
+| Qwen3-TTS cancel latency (core `tts` pack, CPU, 0006), median / p90 / max | 16 / 40 / 49 ms against 15 / 32 / 48 ms |
 
 The first version of 0006 accumulated every tile in F32: F16 GEMM +125–163%
 time, F16-weight pp512 at 0.44 of `v0.8.0-3` and a median TTS cancel latency
 of 65 ms. Blocks of 16 or 32 products per lane, instead of 64, gave pp512
-at 0.79 or 0.87. `linux_armv8.2_2` runs the conversion flush; it is not yet
-measured on that hardware.
+at 0.79 or 0.87.
 
 ## Checks
 
@@ -105,10 +135,32 @@ measured on that hardware.
   that reaches the next row. It reports every failing check. The macOS
   release lanes also run it built for `armv8.2-a+fp16+dotprod` to cover the
   conversion branch.
-- Without the patches it fails on Apple silicon: 16 checks with 0005 alone
-  (GEMM overflow and block cases, eight-query ordinary attention); without
-  0006's flush, recompute or tail block, 1 to 202 checks.
-  A `n - n % 4` tail mutation of 0005 fails 15 checks through the padding. Zero-length and
+- For 0005, on a module with FP16 vector arithmetic and no SVE, the test
+  models the patched single-column dot (upstream's F16 multiply-adds,
+  reduction order and 64-step block, then F32 lanes and the double tail) and
+  requires a bit-exact match wherever the model is finite; elsewhere it
+  requires the exact sum. Its 0005 cases: lane sums exact in F16 within one
+  block but not with a 65th step (a longer block), a product of 2048 that
+  absorbs later products of 1 (F16 rounding, reduction order, a shorter
+  block), one overflowing lane in each half, accumulator and block of one-
+  and multi-block rows,
+  overflows of opposite sign that meet as NaN in the F32 lanes or in the F16
+  reduction (a finiteness check that misses NaN), product overflow at widths
+  around the recompute's eight-element and scalar tails, rows without NaN
+  padding (an over-read the recompute would hide), and NaN and infinite
+  inputs (only there: upstream's other paths assert a finite sum in debug
+  builds).
+- On Apple silicon (FMLAL and conversion builds) it fails 36 checks without
+  0005, 19 with the first 0005 (bit-identity only), and 1 to 40 with each of
+  29 mutants of the new 0005: no flush, an in-line block of 65 or 128 steps,
+  blocks of 32, 63, 65 or 128 steps, a flush that overwrites instead of
+  adding, no recompute, recompute always, finiteness checks of infinity only,
+  the low half or one lane, a dropped scalar tail, partial block or high half,
+  over-reads, a reordered F16 reduction or a dropped accumulator, and a
+  recompute that skips a tail, drops its high half or returns NaN. Moving the
+  recompute's eight-element loop into its scalar tail is equivalent and
+  survives. Without 0006's flush, recompute or tail block it failed 1 to 202
+  checks when 0006 landed; those GEMM cases are unchanged. Zero-length and
   strided query tensors are not covered.
 - `tools/linux_cpu_variant_smoke.c` runs the single-column magnitude sweep on
   each selected Linux arm64 CPU module and rejects nonfinite matmul output.
@@ -128,8 +180,7 @@ is on `main`:
 - Run the test against each published CPU module: Android variants on a
   device, Linux variants, Apple.
 - Compare prefill/decode throughput on an Android device for each variant with
-  the previous release, F16 KV with FA off and on, and an F16-weight model on
-  `linux_armv8.2_2`.
+  the previous release, F16 KV with FA off and on.
 
 Retire a patch when the pinned upstream and every compatibility lane carry a
 verified upstream repair; then add its `upstream_fix` text.

patches/llama.cpp/0005-arm-f16-dot-f32-accumulation.patch

@@ -1,18 +1,14 @@
 --- a/ggml/src/ggml-cpu/vec.cpp
 +++ b/ggml/src/ggml-cpu/vec.cpp
-@@ -263,4 +263,5 @@
+@@ -262,5 +262,93 @@
+ }
 
- void ggml_vec_dot_f16(int n, float * GGML_RESTRICT s, size_t bs, ggml_fp16_t * GGML_RESTRICT x, size_t bx, ggml_fp16_t * GGML_RESTRICT y, size_t by, int nrc) {
-+    GGML_ASSERT(n >= 0 && "llamadart F16 dot F32 accumulation");
-     assert(nrc == 1);
-     GGML_UNUSED(nrc);
-@@ -272,5 +273,48 @@
-
-
--#if defined(GGML_SIMD)
 +#if defined(__ARM_NEON) && defined(__ARM_FEATURE_FMA) && defined(__ARM_FP16_FORMAT_IEEE) && \
 +    defined(__ARM_FEATURE_FP16_VECTOR_ARITHMETIC) && !defined(__ARM_FEATURE_SVE)
-+    // F16 products and partial sums can overflow before the F32 result is stored.
++#define LLAMADART_F16_DOT_BLOCKED
++// F32 products and sums, for a dot product whose F16 sums were nonfinite.
++__attribute__((noinline))
++static void ggml_vec_dot_f16_f32_sums(int n, float * GGML_RESTRICT s, const ggml_fp16_t * GGML_RESTRICT x, const ggml_fp16_t * GGML_RESTRICT y) {
 +    float32x4_t sum_lo[4], sum_hi[4];
 +    for (int j = 0; j < 4; ++j) {
 +        sum_lo[j] = vdupq_n_f32(0.0f);
@@ -49,10 +45,78 @@
 +        sum_lo[0] = vaddq_f32(sum_lo[0], sum_lo[j]);
 +        sum_hi[0] = vaddq_f32(sum_hi[0], sum_hi[j]);
 +    }
-+    sumf = (ggml_float) vaddvq_f32(vaddq_f32(sum_lo[0], sum_hi[0]));
++    ggml_float sumf = (ggml_float) vaddvq_f32(vaddq_f32(sum_lo[0], sum_hi[0]));
 +    for (; i < n; ++i) {
 +        sumf += (ggml_float) GGML_CPU_FP16_TO_FP32(x[i]) * (ggml_float) GGML_CPU_FP16_TO_FP32(y[i]);
 +    }
++    *s = sumf;
++}
++
++// Upstream's F16 multiply-adds over elements [i, end) and its F16 reduction.
++static inline float16x8_t ggml_vec_dot_f16_block(const ggml_fp16_t * GGML_RESTRICT x, const ggml_fp16_t * GGML_RESTRICT y, int i, int end) {
++    float16x8_t sum[4] = { vdupq_n_f16(0.0f), vdupq_n_f16(0.0f), vdupq_n_f16(0.0f), vdupq_n_f16(0.0f) };
++    for (; i < end; i += 32) {
++        for (int j = 0; j < 4; ++j) {
++            sum[j] = vfmaq_f16(sum[j], vld1q_f16((const __fp16 *) (x + i + 8*j)),
++                                       vld1q_f16((const __fp16 *) (y + i + 8*j)));
++        }
++    }
++    sum[0] = vaddq_f16(sum[0], sum[2]);
++    sum[1] = vaddq_f16(sum[1], sum[3]);
++    return vaddq_f16(sum[0], sum[1]);
++}
++
++// Rows longer than one block add each block's F16 sums into F32.
++__attribute__((noinline))
++static void ggml_vec_dot_f16_blocks(int n, float * GGML_RESTRICT s, const ggml_fp16_t * GGML_RESTRICT x, const ggml_fp16_t * GGML_RESTRICT y) {
++    const int np = n & ~31;
++    float32x4_t sum_lo = vdupq_n_f32(0.0f);
++    float32x4_t sum_hi = vdupq_n_f32(0.0f);
++    for (int i = 0; i < np; i += 64 * 32) {
++        const float16x8_t sum = ggml_vec_dot_f16_block(x, y, i, np - i < 64 * 32 ? np : i + 64 * 32);
++        sum_lo = vaddq_f32(sum_lo, vcvt_f32_f16(vget_low_f16(sum)));
++        sum_hi = vaddq_f32(sum_hi, vcvt_f32_f16(vget_high_f16(sum)));
++    }
++    const float sum = vaddvq_f32(vaddq_f32(sum_lo, sum_hi));
++    if (__builtin_expect(!isfinite(sum), 0)) {
++        ggml_vec_dot_f16_f32_sums(n, s, x, y);
++        return;
++    }
++    ggml_float sumf = sum;
++    for (int i = np; i < n; ++i) {
++        sumf += (ggml_float)(GGML_CPU_FP16_TO_FP32(x[i])*GGML_CPU_FP16_TO_FP32(y[i]));
++    }
++    *s = sumf;
++}
++#endif
++
+ void ggml_vec_dot_f16(int n, float * GGML_RESTRICT s, size_t bs, ggml_fp16_t * GGML_RESTRICT x, size_t bx, ggml_fp16_t * GGML_RESTRICT y, size_t by, int nrc) {
++    GGML_ASSERT(n >= 0 && "llamadart F16 dot F32 accumulation");
+     assert(nrc == 1);
+     GGML_UNUSED(nrc);
+@@ -272,5 +360,24 @@
+
+
+-#if defined(GGML_SIMD)
++#if defined(LLAMADART_F16_DOT_BLOCKED)
++    // F16 sums overflow past 65504. Upstream's F16 multiply-adds and reduction
++    // run in its order on at most 64 products per lane, then add in F32; a
++    // nonfinite F16 or F32 sum is recomputed with F32 products and sums.
++    const int np = n & ~31;
++    if (np > 64 * 32) {
++        ggml_vec_dot_f16_blocks(n, s, x, y);
++        return;
++    }
++    const float16x8_t sum = ggml_vec_dot_f16_block(x, y, 0, np);
++    // sum - sum is zero in every lane only when every lane is finite.
++    if (__builtin_expect(vminvq_u16(vceqzq_f16(vsubq_f16(sum, sum))) == 0, 0)) {
++        ggml_vec_dot_f16_f32_sums(n, s, x, y);
++        return;
++    }
++    sumf = (ggml_float) vaddvq_f32(vaddq_f32(vcvt_f32_f16(vget_low_f16(sum)), vcvt_f32_f16(vget_high_f16(sum))));
++    for (int i = np; i < n; ++i) {
++        sumf += (ggml_float)(GGML_CPU_FP16_TO_FP32(x[i])*GGML_CPU_FP16_TO_FP32(y[i]));
++    }
 +#elif defined(GGML_SIMD)
      #if defined(__ARM_FEATURE_SVE)
          const int ggml_f16_epr = svcnth();

patches/llama.cpp/README.md

@@ -26,10 +26,10 @@ arithmetic or llamafile is a build workaround that slows unrelated kernels
 than a few lines.
 
 Published ARM64 artifacts return nonfinite results for finite F16 dot
-products. Patch 0005 changes only the non-SVE NEON `ggml_vec_dot_f16`
-(FMLAL/FMLAL2 with FP16 FML, conversion and F32 FMA otherwise); patch 0006
-makes the FP16 llamafile GEMM tile add its F16 sums into F32 every 512 K
-elements and redo in F32 any tile with a nonfinite result. Neither changes ISA
+products. Patch 0005 makes the non-SVE NEON `ggml_vec_dot_f16` add its F16
+sums into F32 every 2048 elements (64 products per lane) and recompute in F32
+any result that is nonfinite; patch 0006 does the same for the FP16 llamafile
+GEMM tile every 512 K elements. Neither changes ISA
 flags, ABI, variants or cache defaults. Performance cost, checks and
 publication gates are in
 [native174_f16_dot_accumulation.md](../../docs/native174_f16_dot_accumulation.md).