diff --git a/docs/linux_arm64_cpu_variants.md b/docs/linux_arm64_cpu_variants.md index feb189d..04af198 100644 --- a/docs/linux_arm64_cpu_variants.md +++ b/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. diff --git a/docs/native174_f16_dot_accumulation.md b/docs/native174_f16_dot_accumulation.md index c217df2..c26c658 100644 --- a/docs/native174_f16_dot_accumulation.md +++ b/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. diff --git a/patches/llama.cpp/0005-arm-f16-dot-f32-accumulation.patch b/patches/llama.cpp/0005-arm-f16-dot-f32-accumulation.patch index 8961154..9a5dfe8 100644 --- a/patches/llama.cpp/0005-arm-f16-dot-f32-accumulation.patch +++ b/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(); diff --git a/patches/llama.cpp/README.md b/patches/llama.cpp/README.md index 239fedb..dacb43c 100644 --- a/patches/llama.cpp/README.md +++ b/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). diff --git a/tests/f16_dot_accumulation_test.cpp b/tests/f16_dot_accumulation_test.cpp index 88b9604..03f1c16 100644 --- a/tests/f16_dot_accumulation_test.cpp +++ b/tests/f16_dot_accumulation_test.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include namespace { @@ -26,6 +27,61 @@ float rounded_f16(float value) { return ggml_fp16_to_fp32(ggml_fp32_to_fp16(value)); } +// Rounds to the nearest F16, ties to even, overflowing to infinity. Products +// and sums of the F16 test inputs below are exact in double, so this models one +// F16 multiply-add or add. +double round_f16(double value) { + if (!std::isfinite(value) || value == 0.0) { + return value; + } + int exponent; + std::frexp(value, &exponent); + const int quantum = std::max(exponent - 11, -24); + const double rounded = std::ldexp(std::nearbyint(std::ldexp(value, -quantum)), quantum); + return std::fabs(rounded) > 65504.0 ? std::copysign(INFINITY, value) : rounded; +} + +// Whether the CPU module multiplies and adds F16 dot products in F16: FP16 +// vector arithmetic without SVE (upstream widens to F32 under SVE). +bool f16_sums = false; + +// Patch 0005 keeps upstream's single-column F16 dot where it is finite: four +// accumulators of eight lanes take one multiply-add per 32-element step, for +// at most 64 steps; then accumulators 0+2 and 1+3 and the two results add in +// F16, and each lane adds into F32. The lanes reduce as (0+4 + 1+5) + (2+6 + +// 3+7) in F32 and elements past the last full step add in double. A nonfinite +// result (returned as such) is recomputed in F32 by the kernel. +constexpr int dot_step = 32; +constexpr int dot_block_steps = 64; +constexpr int dot_block = dot_step * dot_block_steps; + +double blocked_f16_dot(const std::vector &k, const std::vector &q) { + const int n = static_cast(k.size()); + const int np = n - n % dot_step; + float lanes[8] = {}; + for (int start = 0; start < np; start += dot_block) { + const int end = std::min(np, start + dot_block); + double sums[dot_step] = {}; + for (int i = start; i < end; ++i) { + sums[i % dot_step] = round_f16(sums[i % dot_step] + k[i] * q[i]); + } + for (int lane = 0; lane < 8; ++lane) { + const double sum02 = round_f16(sums[lane] + sums[16 + lane]); + const double sum13 = round_f16(sums[8 + lane] + sums[24 + lane]); + lanes[lane] += static_cast(round_f16(sum02 + sum13)); + } + } + float pairs[4]; + for (int lane = 0; lane < 4; ++lane) { + pairs[lane] = lanes[lane] + lanes[4 + lane]; + } + double sum = (pairs[0] + pairs[1]) + (pairs[2] + pairs[3]); + for (int i = np; i < n; ++i) { + sum += k[i] * q[i]; + } + return sum; +} + class Graph { public: explicit Graph(ggml_backend_t backend) : backend_(backend) { @@ -63,8 +119,14 @@ class Graph { bool check_result(float actual, double expected, double tolerance, const char *name, int width, int cols) { const double error = std::fabs(static_cast(actual) - expected); - if (!std::isfinite(actual) || !std::isfinite(expected) || - !std::isfinite(error) || error > tolerance) { + bool pass = std::isfinite(actual) && std::isfinite(expected) && + std::isfinite(error) && error <= tolerance; + if (std::isnan(expected)) { + pass = std::isnan(actual); + } else if (std::isinf(expected)) { + pass = actual == expected; + } + if (!pass) { std::fprintf(stderr, "FAIL: %s width=%d cols=%d expected=%.17g actual=%.9g tolerance=%.9g\n", name, width, cols, expected, actual, tolerance); @@ -74,29 +136,18 @@ bool check_result(float actual, double expected, double tolerance, return true; } -// Patch 0006 sums each F16 GEMM lane over at most this many K elements (eight -// lanes, so 64 products per lane) before adding it in F32. -constexpr int gemm_block = 512; -constexpr int gemm_block_steps = gemm_block / 8; - -enum class Pattern { - positive, cancellation, large_product, maximum, subnormal, - // Every lane of every GEMM block sums exactly to 65504, the F16 maximum. - block_maximum, - // As block_maximum but each lane sums to 65536, which overflows in F16. - block_overflow, - // Each lane sums 64 products to 2047, exact in F16; a 65th product of 32 - // gives the odd 2079, which F16 cannot hold. - block_exact, -}; - // Sixteen identical weight rows times `cols` identical query columns. One // column takes the vector dot kernel; eight or more take llamafile GEMM where -// it is built. NaN padding after each weight row fails any read past `width`. -void dot_case(ggml_backend_t backend, int width, float magnitude, - Pattern pattern, ggml_type weight_type, const char *name, int cols = 1) { +// it is built. NaN padding after each weight row fails any read past the row, +// unless a nonfinite result is recomputed: without padding such a read takes +// the next row's finite values instead. One F16 column must match +// blocked_f16_dot exactly where that is finite; everything else must match +// the exact sum. +void values_case(ggml_backend_t backend, const std::vector &key_values, + const std::vector &query_values, ggml_type weight_type, + const char *name, int cols = 1, int padding = 8) { constexpr int rows = 16; - constexpr int padding = 8; + const int width = static_cast(key_values.size()); const int stride = width + padding; Graph graph(backend); auto *padded = ggml_new_tensor_2d(graph.ctx, weight_type, stride, rows); @@ -109,30 +160,12 @@ void dot_case(ggml_backend_t backend, int width, float magnitude, std::vector keys(static_cast(stride) * rows, NAN); std::vector half_keys(keys.size(), ggml_fp32_to_fp16(NAN)); std::vector queries(static_cast(width) * cols); + std::vector actual_k(width), actual_q(width); double reference = 0.0; double absolute_sum = 0.0; for (int i = 0; i < width; ++i) { - float kv = magnitude * (1.0f + 0.1f * (i % 5)); - float qv = magnitude * (1.0f + 0.1f * ((i * 3) % 5)); - if (pattern == Pattern::cancellation) { - kv = magnitude; - qv = i < width / 2 ? magnitude : -magnitude; - } else if (pattern == Pattern::large_product) { - kv = qv = 256.0f; - } else if (pattern == Pattern::maximum) { - kv = qv = 65504.0f; - } else if (pattern == Pattern::subnormal) { - kv = ggml_fp16_to_fp32(static_cast(1 + i % 7)); - qv = i % 2 == 0 ? 1.0f : -1.0f; - } else if (pattern == Pattern::block_maximum || pattern == Pattern::block_overflow) { - const bool last_step = (i % gemm_block) / 8 == gemm_block_steps - 1; - kv = 32.0f; - qv = 2048.0f / gemm_block_steps - - (pattern == Pattern::block_maximum && last_step ? 1.0f : 0.0f); - } else if (pattern == Pattern::block_exact) { - kv = (i % gemm_block) / 8 == 1 ? 31.0f : 32.0f; - qv = 1.0f; - } + const float kv = key_values[i]; + const float qv = query_values[i]; for (int r = 0; r < rows; ++r) { keys[static_cast(r) * stride + i] = kv; half_keys[static_cast(r) * stride + i] = ggml_fp32_to_fp16(kv); @@ -140,13 +173,22 @@ void dot_case(ggml_backend_t backend, int width, float magnitude, for (int c = 0; c < cols; ++c) { queries[static_cast(c) * width + i] = qv; } - const double actual_k = weight_type == GGML_TYPE_F16 ? rounded_f16(kv) : kv; + actual_k[i] = weight_type == GGML_TYPE_F16 ? rounded_f16(kv) : kv; // The F16 CPU type trait rounds the F32 query before its dot product. - const double actual_q = weight_type == GGML_TYPE_F16 ? rounded_f16(qv) : qv; - const double product = actual_k * actual_q; + actual_q[i] = weight_type == GGML_TYPE_F16 ? rounded_f16(qv) : qv; + const double product = actual_k[i] * actual_q[i]; reference += product; absolute_sum += std::fabs(product); } + double expected = reference; + double tolerance = 1e-7 + 2e-6 * absolute_sum; + if (cols == 1 && weight_type == GGML_TYPE_F16 && f16_sums) { + const double blocked = blocked_f16_dot(actual_k, actual_q); + if (std::isfinite(blocked)) { + expected = static_cast(blocked); + tolerance = 0.0; + } + } if (weight_type == GGML_TYPE_F16) { ggml_backend_tensor_set(padded, half_keys.data(), 0, half_keys.size() * sizeof(ggml_fp16_t)); } else { @@ -157,12 +199,93 @@ void dot_case(ggml_backend_t backend, int width, float magnitude, std::vector actual(static_cast(rows) * cols); ggml_backend_tensor_get(out, actual.data(), 0, actual.size() * sizeof(float)); for (float value : actual) { - if (!check_result(value, reference, 1e-7 + 2e-6 * absolute_sum, name, width, cols)) { + if (!check_result(value, expected, tolerance, name, width, cols)) { break; } } } +// Patch 0006 sums each F16 GEMM lane over at most this many K elements (eight +// lanes, so 64 products per lane) before adding it in F32. +constexpr int gemm_block = 512; +constexpr int gemm_block_steps = gemm_block / 8; + +enum class Pattern { + positive, cancellation, large_product, maximum, subnormal, + // Every lane of every GEMM block sums exactly to 65504, the F16 maximum. + block_maximum, + // As block_maximum but each lane sums to 65536, which overflows in F16. + block_overflow, + // Each lane sums 64 products to 2047, exact in F16; a 65th product of 32 + // gives the odd 2079, which F16 cannot hold. + block_exact, + // block_exact for the single-column dot's 64 steps of 32 elements. + dot_block_exact, + // Lanes 0-7 of the first step take 2048, after which each product of 1 + // rounds away in F16; with the other accumulators' lanes this pins + // upstream's F16 multiply-add, reduction order and block length. + f16_rounding, +}; + +void dot_case(ggml_backend_t backend, int width, float magnitude, Pattern pattern, + ggml_type weight_type, const char *name, int cols = 1, int padding = 8) { + std::vector keys(width), queries(width); + for (int i = 0; i < width; ++i) { + float kv = magnitude * (1.0f + 0.1f * (i % 5)); + float qv = magnitude * (1.0f + 0.1f * ((i * 3) % 5)); + if (pattern == Pattern::cancellation) { + kv = magnitude; + qv = i < width / 2 ? magnitude : -magnitude; + } else if (pattern == Pattern::large_product) { + kv = qv = 256.0f; + } else if (pattern == Pattern::maximum) { + kv = qv = 65504.0f; + } else if (pattern == Pattern::subnormal) { + kv = ggml_fp16_to_fp32(static_cast(1 + i % 7)); + qv = i % 2 == 0 ? 1.0f : -1.0f; + } else if (pattern == Pattern::block_maximum || pattern == Pattern::block_overflow) { + const bool last_step = (i % gemm_block) / 8 == gemm_block_steps - 1; + kv = 32.0f; + qv = 2048.0f / gemm_block_steps - + (pattern == Pattern::block_maximum && last_step ? 1.0f : 0.0f); + } else if (pattern == Pattern::block_exact) { + kv = (i % gemm_block) / 8 == 1 ? 31.0f : 32.0f; + qv = 1.0f; + } else if (pattern == Pattern::dot_block_exact) { + kv = (i % dot_block) / dot_step == 1 ? 31.0f : 32.0f; + qv = 1.0f; + } else if (pattern == Pattern::f16_rounding) { + kv = i < 8 ? 64.0f : 1.0f; + qv = i < 8 ? 32.0f : 1.0f; + } + keys[i] = kv; + queries[i] = qv; + } + values_case(backend, keys, queries, weight_type, name, cols, padding); +} + +// One F16 multiply-add of 256 * 256 = 65536 overflows the lane at `index` +// (step index / 32, position index % 32); a second index takes -65536. All +// other products are 1, so every exact result is finite. +void lane_overflow_case(ggml_backend_t backend, int width, int index, const char *name, + int negative_index = -1) { + std::vector keys(width, 1.0f), queries(width, 1.0f); + keys[index] = queries[index] = 256.0f; + if (negative_index >= 0) { + keys[negative_index] = 256.0f; + queries[negative_index] = -256.0f; + } + values_case(backend, keys, queries, GGML_TYPE_F16, name); +} + +// A NaN or infinite input stays nonfinite, in the F16 steps and in the tail. +void nonfinite_input_case(ggml_backend_t backend, int width, int index, float value, + const char *name) { + std::vector keys(width, 1.0f), queries(width, 1.0f); + keys[index] = value; + values_case(backend, keys, queries, GGML_TYPE_F16, name); +} + // Sixteen weight rows times twelve query columns with exact integer products. // In the K tail after one GEMM block, lane 0 overflows F16 for rows 5-6 times // columns 8-9 and lane 1 for rows 1 and 7 times columns 2 and 11. The 4x6 tiles @@ -320,10 +443,16 @@ int main() { ggml_backend_reg_get_proc_address(reg, "ggml_backend_get_features")); require(get_features != nullptr, "CPU features unavailable"); std::printf("CPU features:"); + bool fp16_va = false; + bool sve = false; for (auto *feature = get_features(reg); feature->name != nullptr; ++feature) { std::printf(" %s=%s", feature->name, feature->value); + fp16_va |= std::strcmp(feature->name, "FP16_VA") == 0; + sve |= std::strcmp(feature->name, "SVE") == 0; } - std::printf("\n"); + f16_sums = fp16_va && !sve; + std::printf("\nSingle-column F16 dot expectation: %s\n", + f16_sums ? "upstream F16 sums in 64-step blocks" : "exact"); for (int magnitude = 30; magnitude <= 100; magnitude += 10) { dot_case(backend, 64, static_cast(magnitude), Pattern::positive, @@ -340,6 +469,53 @@ int main() { dot_case(backend, 32, 0.0f, Pattern::large_product, GGML_TYPE_F16, "large-product"); dot_case(backend, 64, 0.0f, Pattern::maximum, GGML_TYPE_F16, "maximum-finite"); dot_case(backend, 64, 100.0f, Pattern::positive, GGML_TYPE_F32, "f32-control"); + for (int width : {33, 39, 40, 41, 63, 65, 127, 129}) { + dot_case(backend, width, 0.0f, Pattern::maximum, GGML_TYPE_F16, "maximum-tail"); + dot_case(backend, width, 1.0f, Pattern::positive, GGML_TYPE_F16, "unpadded-tail", 1, 0); + } + for (int width : {dot_block - dot_step, dot_block, dot_block + dot_step, dot_block + 39, + 2 * dot_block + dot_step + 7}) { + dot_case(backend, width, 0.0f, Pattern::dot_block_exact, GGML_TYPE_F16, "dot-block-exact"); + } + for (int width : {dot_block - 1, dot_block + dot_step - 1, 2 * dot_block + 5}) { + dot_case(backend, width, 0.0f, Pattern::f16_rounding, GGML_TYPE_F16, "f16-rounding"); + } + { + // Rows of one block: acc0 lane 7 (high half) of the first step and acc3 + // lane 7 of the second; overflows of opposite sign in acc0 lanes 0 and 4, + // and in acc0 and acc2 lane 0. + constexpr int width = 2 * dot_step + 5; + for (int index : {7, dot_step + 31}) { + lane_overflow_case(backend, width, index, "short-lane-overflow"); + } + lane_overflow_case(backend, width, 0, "short-split-overflow", dot_step + 4); + lane_overflow_case(backend, width, 0, "short-split-overflow", dot_step + 16); + } + { + constexpr int width = 2 * dot_block + dot_step + 5; + // The first block's acc0 lane 0, the second's acc0 lane 7 (high half), + // the last step of the first block in acc2 lane 4, and the final partial + // block's acc3 lane 7. + for (int index : {3 * dot_step, (dot_block_steps + 5) * dot_step + 7, + (dot_block_steps - 1) * dot_step + 20, 2 * dot_block + 31}) { + lane_overflow_case(backend, width, index, "lane-overflow"); + } + // +inf in lane 0 and -inf in lane 4 meet in the F32 sum; +inf in acc0 and + // -inf in acc2 meet in the F16 reduction. Both give NaN, not infinity. + lane_overflow_case(backend, width, 3 * dot_step, "split-overflow-f32", + (dot_block_steps + 6) * dot_step + 4); + lane_overflow_case(backend, width, 3 * dot_step, "split-overflow-f16", + 4 * dot_step + 16); + } + // Upstream's other F16 dot paths assert a finite sum in debug builds. + if (f16_sums) { + for (int width : {2 * dot_step + 5, 2 * dot_block + dot_step + 5}) { + for (int index : {9, width - 3}) { + nonfinite_input_case(backend, width, index, NAN, "nan-input"); + nonfinite_input_case(backend, width, index, INFINITY, "inf-input"); + } + } + } dot_case(backend, 256, 128.0f, Pattern::cancellation, GGML_TYPE_F16, "hidden-overflow"); for (int cols : {8, 16, 32}) { dot_case(backend, 64, 100.0f, Pattern::positive, GGML_TYPE_F16, "gemm-positive", cols); diff --git a/tools/linux_cpu_variant_smoke.c b/tools/linux_cpu_variant_smoke.c index e52c232..cb72a93 100644 --- a/tools/linux_cpu_variant_smoke.c +++ b/tools/linux_cpu_variant_smoke.c @@ -170,7 +170,9 @@ static int run_finite_dot_cases(ggml_backend_dev_t device) { float actual; ggml_backend_tensor_get(out, &actual, 0, sizeof(actual)); const double error = fabs((double)actual - reference); - if (!isfinite(actual) || !isfinite(error) || error > 1e-7 + 2e-6 * absolute_sum) { + // Patch 0005 keeps upstream's F16 multiply-adds and reduction: at most + // four F16 roundings, each within 2^-11 of the absolute sum, for 64 elements. + if (!isfinite(actual) || !isfinite(error) || error > 1e-7 + 0x1p-9 * absolute_sum) { fprintf(stderr, "finite dot case %d mag %.0f expected %.17g got %.9g\n", test_case, magnitude, reference, actual); status = 1;