Pin F16 dot negative overflow and bit-identity in the 0005 test

docs/native174_f16_dot_accumulation.md

@@ -139,11 +139,16 @@ at 0.79 or 0.87.
   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
+  requires a bit-exact match to the first 0005's F32 recompute (fused
+  multiply-adds in four accumulators, then the accumulators, halves and pairs
+  added in order, and a double tail). 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,
+  and multi-block rows, the same overflow in the negative direction (a check
+  that sees only positive overflow), mixed-exponent non-integer rows whose F32
+  rounding makes the order of the sums visible, finite (F32 lane order, a block
+  tail summed in float) or recomputed (accumulator order),
   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

tests/f16_dot_accumulation_test.cpp

@@ -82,6 +82,37 @@ double blocked_f16_dot(const std::vector<double> &k, const std::vector<double> &
   return sum;
 }
 
+// The F32 recompute of a nonfinite row, in the order of the first version of
+// 0005: fused F32 multiply-adds in four accumulators of eight lanes per
+// 32-element step, then single steps of eight into accumulator 0, the
+// accumulators added in order, the low and high four lanes added, the four
+// results added pairwise, and a double tail.
+float f32_recompute_dot(const std::vector<double> &k, const std::vector<double> &q) {
+  const int n = static_cast<int>(k.size());
+  const int np = n - n % dot_step;
+  const int np8 = n - n % 8;
+  float acc[4][8] = {};
+  int i = 0;
+  for (; i < np8; ++i) {
+    float &lane = i < np ? acc[(i % dot_step) / 8][i % 8] : acc[0][i % 8];
+    lane = std::fmaf(static_cast<float>(k[i]), static_cast<float>(q[i]), lane);
+  }
+  for (int j = 1; j < 4; ++j) {
+    for (int lane = 0; lane < 8; ++lane) {
+      acc[0][lane] += acc[j][lane];
+    }
+  }
+  float halves[4];
+  for (int lane = 0; lane < 4; ++lane) {
+    halves[lane] = acc[0][lane] + acc[0][4 + lane];
+  }
+  double sum = (halves[0] + halves[1]) + (halves[2] + halves[3]);
+  for (; i < n; ++i) {
+    sum += k[i] * q[i];
+  }
+  return static_cast<float>(sum);
+}
+
 class Graph {
  public:
   explicit Graph(ggml_backend_t backend) : backend_(backend) {
@@ -141,8 +172,8 @@ bool check_result(float actual, double expected, double tolerance,
 // 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.
+// blocked_f16_dot exactly where that is finite and f32_recompute_dot exactly
+// where it is not; everything else must match the exact sum.
 void values_case(ggml_backend_t backend, const std::vector<float> &key_values,
                  const std::vector<float> &query_values, ggml_type weight_type,
                  const char *name, int cols = 1, int padding = 8) {
@@ -184,10 +215,9 @@ void values_case(ggml_backend_t backend, const std::vector<float> &key_values,
   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<float>(blocked);
-      tolerance = 0.0;
-    }
+    expected = std::isfinite(blocked) ? static_cast<float>(blocked)
+                                      : f32_recompute_dot(actual_k, actual_q);
+    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));
@@ -265,19 +295,47 @@ void dot_case(ggml_backend_t backend, int width, float magnitude, Pattern patter
 }
 
 // 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.
+// (step index / 32, position index % 32); a second index takes -65536, and
+// `negative` makes the first -65536 as well. 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) {
+                        int negative_index = -1, bool negative = false) {
   std::vector<float> keys(width, 1.0f), queries(width, 1.0f);
-  keys[index] = queries[index] = 256.0f;
+  keys[index] = 256.0f;
+  queries[index] = negative ? -256.0f : 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);
 }
 
+// Products of mixed binary exponents with fractions that F32 sums round, so a
+// changed order of the F32 additions changes the result. The lane `i % 8`
+// takes the exponent. Each element past the last full step is 3 * 2^-11: under
+// half an ulp of the row's sum (about 10^5 at width 4133), but five of them
+// are not, so a tail summed in float differs from one summed in double.
+// `overflow_index` >= 0 adds a product of 65536 there.
+void spread_case(ggml_backend_t backend, int width, int overflow_index, const char *name) {
+  constexpr int exponents[8] = {5, -7, 3, -9, 6, -6, 2, -10};
+  const int np = width - width % dot_step;
+  std::vector<float> keys(width), queries(width);
+  for (int i = 0; i < width; ++i) {
+    if (i >= np) {
+      keys[i] = 0x3p-8f;
+      queries[i] = 0x1p-3f;
+      continue;
+    }
+    const int exponent = exponents[i % 8];
+    keys[i] = (1.0f + ((7 * i + 3) % 11) / 16.0f) * std::ldexp(1.0f, exponent / 2);
+    queries[i] = (1.0f + ((5 * i + 1) % 13) / 32.0f) * std::ldexp(1.0f, exponent - exponent / 2);
+  }
+  if (overflow_index >= 0) {
+    keys[overflow_index] = queries[overflow_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) {
@@ -490,6 +548,9 @@ int main() {
     }
     lane_overflow_case(backend, width, 0, "short-split-overflow", dot_step + 4);
     lane_overflow_case(backend, width, 0, "short-split-overflow", dot_step + 16);
+    for (int index : {7, dot_step + 31, width - 3}) {
+      lane_overflow_case(backend, width, index, "short-negative-overflow", -1, true);
+    }
   }
   {
     constexpr int width = 2 * dot_block + dot_step + 5;
@@ -506,6 +567,21 @@ int main() {
                        (dot_block_steps + 6) * dot_step + 4);
     lane_overflow_case(backend, width, 3 * dot_step, "split-overflow-f16",
                        4 * dot_step + 16);
+    for (int index : {3 * dot_step, 2 * dot_block + 31}) {
+      lane_overflow_case(backend, width, index, "negative-overflow", -1, true);
+    }
+  }
+  // Elsewhere F32 sums of these rows drift past the exact-sum tolerance.
+  if (f16_sums) {
+    constexpr int short_width = 2 * dot_step + 5;
+    constexpr int long_width = 2 * dot_block + dot_step + 5;
+    spread_case(backend, 64, -1, "spread");
+    spread_case(backend, short_width, -1, "spread-tail");
+    spread_case(backend, short_width, 11, "spread-recompute");
+    spread_case(backend, short_width + 3, dot_step + 20, "spread-recompute-tail");
+    spread_case(backend, long_width, -1, "spread-blocks");
+    spread_case(backend, long_width + 3, -1, "spread-blocks-tail");
+    spread_case(backend, long_width, 5 * dot_step + 3, "spread-blocks-recompute");
   }
   // Upstream's other F16 dot paths assert a finite sum in debug builds.
   if (f16_sums) {

tools/linux_cpu_variant_smoke.c

@@ -171,8 +171,10 @@ static int run_finite_dot_cases(ggml_backend_dev_t device) {
       ggml_backend_tensor_get(out, &actual, 0, sizeof(actual));
       const double error = fabs((double)actual - reference);
       // 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) {
+      // four F16 roundings, each within 2^-11 of the absolute sum, for 64
+      // elements. The F32 control runs no patched code.
+      const double tolerance = 1e-7 + (f32_control ? 2e-6 : 0x1p-9) * absolute_sum;
+      if (!isfinite(actual) || !isfinite(error) || error > tolerance) {
         fprintf(stderr, "finite dot case %d mag %.0f expected %.17g got %.9g\n",
                 test_case, magnitude, reference, actual);
         status = 1;