HMMA with an f16 accumulator rounds once per 8 products, as Ampere does (chain change)

dm/gfx/nvn/dmsh_chain.h

@@ -1,3 +1,3 @@
 // Generated file (dm/gfx/tools/dmsh_chain.py --header): the B shader library translation chain version this source tree expects (default options), one for every game;
 // the host loads only <library dir>/b-<V>/. Run dmsh_chain.py --header after changing the chain's tools and commit this file with the change
-#define DMSH_CHAIN "3df44a5a4e0c"
+#define DMSH_CHAIN "c15175f68cd7"

dm/gfx/tests/test_dmsh_hmma.py

@@ -154,6 +154,7 @@
 #import <Metal/Metal.h>
 static float h2f(uint16_t h) { __fp16 v; memcpy(&v, &h, 2); return (float)v; }
 static uint16_t f2h(float f) { __fp16 v = (__fp16)f; uint16_t h; memcpy(&h, &v, 2); return h; }
+static float r16(float f) { return h2f(f2h(f)); }   // one rounding to f16 (RNE)
 static id<MTLComputePipelineState> pipeline(id<MTLDevice> dev, const char *path, NSString *name) {
     NSError *e = nil; NSData *lib_data = [NSData dataWithContentsOfFile:@(path)];
     id<MTLLibrary> lib = [dev newLibraryWithData:dispatch_data_create(lib_data.bytes, lib_data.length, nil, DISPATCH_DATA_DESTRUCTOR_DEFAULT) error:&e];
@@ -168,8 +169,8 @@
     [ce dispatchThreadgroups:MTLSizeMake(groups, 1, 1) threadsPerThreadgroup:MTLSizeMake(128, 1, 1)]; [ce endEncoding]; [cb commit]; [cb waitUntilCompleted];
 }
 '''
-# The GPU check: 64 threadgroups of 128 lanes, random halves in [-1, 1]; the reference sums in float and rounds once; the matrix hardware rounds the f16 accumulator as
-# it goes, as Ampere's f16 accumulator does, so the tolerance is 0.01 + 0.4 % of the value
+# The GPU check: 64 threadgroups of 128 lanes, random halves in [-1, 1]; the reference is Ampere's: C, then k 0-7 summed in float in order and rounded once to f16,
+# then k 8-15 the same way (dmsh_build.py HMMA_BODY), and every output must match it bit for bit
 RUNNER = HALF + r'''
 int main(int argc, char **argv) {
     @autoreleasepool {
@@ -194,19 +195,19 @@
                 int g = L / 4, t = L % 4; const uint32_t *r = o + (warp * 32 + L) * 2;
                 for (int q2 = 0; q2 < 2; q2++) for (int j = 0; j < 2; j++) {
                     int row = g + 8 * q2, col = 2 * t + j; float ref = C[row][col];
-                    for (int k = 0; k < 16; k++) ref += A[row][k] * B[k][col];
+                    for (int k = 0; k < 16; k++) { ref = fmaf(A[row][k], B[k][col], ref); if (k == 7 || k == 15) ref = r16(ref); }
                     double d = fabs(h2f((r[q2] >> (16 * j)) & 0xffff) - ref);
                     if (d > worst) worst = d;
-                    bad += d > 0.01 + fabs(ref) * 0.004;
+                    bad += (r[q2] >> (16 * j) & 0xffff) != f2h(ref);
                 }
             }
         }
-        printf("%d SIMD-groups, max abs error %.4g, %ld outputs over tolerance\n", n / 32, worst, bad);
+        printf("%d SIMD-groups, max abs error %.4g, %ld outputs that differ from Ampere's\n", n / 32, worst, bad);
         return bad != 0;
     }
 }
 '''
-# The run of calls: E0 = A0 (B0 + B1) + C and E1 = A1 (B2 + B3), each two HMMAs with an f16 accumulator in between: the tolerance is twice the single call's
+# The run of calls: E0 = A0 B0 + A0 B1 + C and E1 = A1 B2 + A1 B3, each two HMMAs with an f16 accumulator in between, bit for bit against the same reference
 SEQ_RUNNER = HALF + r'''
 int main(int argc, char **argv) {
     @autoreleasepool {
@@ -227,14 +228,16 @@
                     int row = g + 8 * q2, col = 2 * t + j;
                     const uint16_t *A = x + 256 * e, *Ba = x + 512 + 256 * e, *Bb = Ba + 128;
                     float ref = e == 0 ? h2f(x[1024 + row * 8 + col]) : 0;
-                    for (int k = 0; k < 16; k++) ref += h2f(A[row * 16 + k]) * (h2f(Ba[k * 8 + col]) + h2f(Bb[k * 8 + col]));
+                    for (int c = 0; c < 2; c++) for (int k = 0; k < 16; k++) {   // two calls, B0 then B1: their k 0-7 and k 8-15 each rounded once
+                        ref = fmaf(h2f(A[row * 16 + k]), h2f((c ? Bb : Ba)[k * 8 + col]), ref); if (k == 7 || k == 15) ref = r16(ref);
+                    }
                     double d = fabs(h2f((r[e * 2 + q2] >> (16 * j)) & 0xffff) - ref);
                     if (d > worst) worst = d;
-                    bad += d > 0.02 + fabs(ref) * 0.008;
+                    bad += (r[e * 2 + q2] >> (16 * j) & 0xffff) != f2h(ref);
                 }
             }
         }
-        printf("%d SIMD-groups, max abs error %.4g, %ld outputs over tolerance\n", warps, worst, bad);
+        printf("%d SIMD-groups, max abs error %.4g, %ld outputs that differ from Ampere's\n", warps, worst, bad);
         return bad != 0;
     }
 }

dm/gfx/tools/dmsh_build.py

@@ -100,15 +100,22 @@ def mul_extended(msl):
     // C and D in the call's own matrices: a slot is read only for a chain's C (cslot) and written only for its D (dslot). A call outside every chain has slot 0 in its
     // flags (no slot bits); working in acc[0..1] it overwrote the accumulator a chain held there between its two calls (Ocarina of Time's agl_ndsr_vr networks: 2,515 such
     // calls in 30 shaders, dotted rows across fog and light, GitHub issue #1537). The flags are constants at each inlined call, so the copies cost nothing
+    // Ampere adds each 8 products to the accumulator and rounds once to f16, k 0-7 then k 8-15 (N_FMA 8 in Khattak and Mikaitis's A100 model, which matches recorded
+    // A100 output in 5,000 of 5,000 cases), where simdgroup_half8x8 rounds after every product (2,420 of them). So each block runs in f32 (each product of two halves
+    // exact, the sum in order) and its result is rounded to f16 (RNE): an f16 C, the first block's result and D are the only roundings
+    // simdgroup_multiply_accumulate takes a float result from half operands and a half C. A lane's thread_elements() are the same two places in a half and a
+    // float matrix, so the rounding is two conversions a lane (a multiply-accumulate with zero A and B into a half result gives the same bits, at twice the cost
+    // of the whole network: NDSR's fa3033d05c0a41d6 0.59 ms against 0.32 on an M5 Max)
     simdgroup_half8x8 b0, b1, c0, c1;
+    simdgroup_float8x8 f0, f1;
     simdgroup_load(b0, h + 256, 16, ulong2(0, 0), true); simdgroup_load(b1, h + 256, 16, ulong2(8, 0), true);
-    if (czero) { simdgroup_multiply(c0, am[0], b0); simdgroup_multiply(c1, am[2], b0); }
-    else {
-        if (cslot) { c0 = acc[k]; c1 = acc[k + 1]; }
-        else { simdgroup_load(c0, h + 384, 8, ulong2(0, 0)); simdgroup_load(c1, h + 384, 8, ulong2(0, 8)); }
-        simdgroup_multiply_accumulate(c0, am[0], b0, c0); simdgroup_multiply_accumulate(c1, am[2], b0, c1);
-    }
-    simdgroup_multiply_accumulate(c0, am[1], b1, c0); simdgroup_multiply_accumulate(c1, am[3], b1, c1);
+    if (czero) { c0 = make_filled_simdgroup_matrix<half, 8, 8>(0.0h); c1 = c0; }
+    else if (cslot) { c0 = acc[k]; c1 = acc[k + 1]; }
+    else { simdgroup_load(c0, h + 384, 8, ulong2(0, 0)); simdgroup_load(c1, h + 384, 8, ulong2(0, 8)); }
+    simdgroup_multiply_accumulate(f0, am[0], b0, c0); simdgroup_multiply_accumulate(f1, am[2], b0, c1);
+    for (int i = 0; i < 2; i++) { c0.thread_elements()[i] = half(f0.thread_elements()[i]); c1.thread_elements()[i] = half(f1.thread_elements()[i]); }   // the first block's rounding
+    simdgroup_multiply_accumulate(f0, am[1], b1, c0); simdgroup_multiply_accumulate(f1, am[3], b1, c1);
+    for (int i = 0; i < 2; i++) { c0.thread_elements()[i] = half(f0.thread_elements()[i]); c1.thread_elements()[i] = half(f1.thread_elements()[i]); }
     if (dslot) { acc[k] = c0; acc[k + 1] = c1; simdgroup_barrier(mem_flags::mem_threadgroup); %(d1)s = 0; return 0; }
     simdgroup_store(c0, h + 384, 8, ulong2(0, 0)); simdgroup_store(c1, h + 384, 8, ulong2(0, 8));
     simdgroup_barrier(mem_flags::mem_threadgroup);

docs/agents/symptoms.md

@@ -960,6 +960,7 @@ One line per symptom in the [subsystem references](../reference/): what you see
 - a Switch 2 game's movie or upscaler network runs slow on the Mac and iPhone, and the GPU profile (`DM_NVN_GPUPROF` with `DM_NVN_DISPATCH_OWN=1`) shows its compute dispatches, small grids of one to four SIMD-groups a threadgroup that loop over tiles, at 1-3 ms each; the dispatches' shaders call `dmsh_hmma16816_f16` hundreds of times (Ampere HMMA through simdgroup_matrix). → [HMMA chains keep their matrices in the SIMD-group](../reference/shaders.md#hmma-chains-keep-their-matrices-in-the-simd-group)
 - a Switch 2 game's movie or upscaler network output has dotted streaks, every other pixel of a few rows wrong, in other places each run, while its input is clean. Through agl's NDSR upscaler a few such pixels a frame grow into white blotches over the whole picture: Ocarina of Time's Kokiri Forest on an M4 Pro on chain 844bf1a5ea8f, with the game's own upscaler (the default); with `DM_UPSCALER=metalfx` in its place, no blotches. → [LDSM reads rows other lanes stored](../reference/shaders.md#ldsm-reads-rows-other-lanes-stored)
 - dotted horizontal rows across the picture at a regular spacing, in fog and light, with the game's upscaler and with `DM_UPSCALER=metalfx` alike (MetalFX replaces only `agl_ndsr.pck`'s shaders): Ocarina of Time in Kokiri Forest, Hyrule Field and the castle town since chain 23a1ea1db092. The rows are the tiles of agl's `agl_ndsr_vr.pck` / `agl_ndsr_vr_light.pck` networks (about 28.8 pixels of the 1536-wide input). → [An HMMA outside a chain overwrote the chain's accumulator](../reference/shaders.md#an-hmma-outside-a-chain-overwrote-the-chains-accumulator)
+- a Switch 2 game's movie or upscaler network gives values a few f16 steps off the console's, in about half the outputs of each HMMA. Nothing looks broken; it is the network's arithmetic, compounded layer by layer. → [HMMA rounds its f16 accumulator once per 8 products](../reference/shaders.md#hmma-rounds-its-f16-accumulator-once-per-8-products)
 - a vertex or compute shader reads zeros from NVN's driver segment `c0`. → [Driver-segment reads outside fragment shaders](../reference/shaders.md#driver-segment-reads-outside-fragment-shaders)
 - flickering hard-edged white tiles (tiled lighting with every tile row getting row 0's or other rows' lights). → [Compute `c0`: dispatch size and uniform buffer table](../reference/shaders.md#compute-c0-dispatch-size-and-uniform-buffer-table)
 - `[RED] shader ... reflection conversion failed: B binding gp_s0 has no resolved slot number`; a geometry shader reading driver segment `c0` reports that its storage-buffer table base was never measured. → [Geometry-stage storage buffer table](../reference/shaders.md#geometry-stage-storage-buffer-table)

docs/reference/shaders.md

@@ -1465,6 +1465,14 @@ Finding, keying and translating shaders (`dm/gfx/tools/shader_split.py`, `dmsh_b
 - **Check:** `test_dmsh_hmma.py` `HmmaOnGpu.test_calls_outside_a_chain_leave_its_slot_alone` (a chain in slot 0 around two calls outside it, on the GPU against a CPU reference): 8,090 outputs over tolerance (largest error 9.6) with the helper before, none after. Ocarina of Time on an M4 Pro, the same Kokiri Forest view: rows on chain 23a1ea1db092 and later; none on 0f510fd5's chain, with the helper's rewrite left out, or with this fix.
 - **Found on:** [The Legend of Zelda: Ocarina of Time](../../dm/games/0400486024DEC000/README.md) (GitHub issue #1537). Counted in the translated MSL (calls outside a chain while slot 0 is held): Ocarina of Time 3,700 in 52 of its 62 HMMA shaders (its NDSR upscaler and movie networks too); [Donkey Kong Bananza](../../dm/games/040074A01BF12000/README.md)'s NDSR 7 in 5 of its 24, so a few of its outputs were wrong as well (not looked at in the game).
 
+### HMMA rounds its f16 accumulator once per 8 products
+
+- **Symptom:** a Switch 2 game's movie or upscaler network gives values a few f16 steps off the console's, in about half the outputs of each HMMA. Nothing looks broken; it is the network's arithmetic, compounded layer by layer.
+- **Cause:** `HMMA_BODY` accumulated on `simdgroup_half8x8`, whose multiply-accumulate rounds the f16 accumulator after every product. Ampere's HMMA with an f16 accumulator adds 8 products to C and rounds once (Khattak and Mikaitis's A100 model, N_FMA 8; [MATLAB-tensor-core](https://github.com/north-numerical-computing/MATLAB-tensor-core)): an HMMA.16816 is k 0-7 rounded to f16, then k 8-15 on top, rounded again. Against the toolbox's 5,000 recorded A100 fp16 cases (C as f16), per-product rounding matches 2,420; summing the 8 products in f32 and rounding once matches all 5,000.
+- **Fix:** `HMMA_BODY` accumulates each block of 8 into a `simdgroup_float8x8` (half operands and a half C, which `simdgroup_multiply_accumulate` takes) and rounds it to f16 per lane through `thread_elements()`. A multiply-accumulate with zero A and B into a half result rounds the same, but took NDSR's network twice as long.
+- **Check:** `test_dmsh_hmma.py`'s GPU runs compare every output bit for bit with the two-block reference (main's body differs in 21,461 of them). Isolated on an M5 Max, NDSR's `fa3033d05c0a41d6` takes 0.316 ms against 0.276 and `d99fbeeb53d38847` 0.446 against 0.393; their outputs change in 1.5 % and 0.8 % of halves.
+- **Found on:** [The Legend of Zelda: Ocarina of Time](../../dm/games/0400486024DEC000/README.md): its NDVE movie network and agl's NDSR upscaler.
+
 ### Driver-segment reads outside fragment shaders
 
 - **Symptom:** a vertex or compute shader reads zeros from NVN's driver segment `c0`.