Merge pull request #40 from itchyfeetleech/perf-exact

engine: SR presets M / L and RR preset F faster, byte-identical

engine/README.md

@@ -1093,6 +1093,23 @@ on the Vulkan shaders unless `[Engine] Network = hip`:
 A model built with no options was compared shader by shader with the explicitly optioned build of each row, and the
 exact switches still reproduce the captures (`build/rr/f/nv/defaults.sh`).
 
+### Byte-identical speed-ups (2026-10-19)
+
+RX 9070 XT, RADV 26.2.2, `d4r-m` 1280x720 -> 2560x1440 with `D4R_ENGINE_M_LAYER_STAMPS=1`. Every output frame, and
+everything `D4R_M_DUMP` writes, is byte-identical to the previous build (M and L, HDR and LDR, the three
+reconstruction tile tiers and the stress inputs).
+
+| change | stage | before -> after (us) |
+|---|---|---|
+| `layer_m.comp`: hardware FP8 saturation where a value cannot be infinite, and the V tiles in a buffer that is dead during the attention (more waves per SIMD) | frame | M 1583 -> 1566, L 2349 -> 2293 |
+| `post_m.comp`: the second history candidate reuses the first one's sample when both have the same vector | reconstruction, uniform motion | 325 -> 288 |
+| `post_m.comp`: both rows of a pixel pair take their Gaussian taps in one pass | reconstruction | 287 -> 265 (L 355 -> 320) |
+| `down_m.comp`: the horizontal pass reads its texels into registers, no shared source tile | downsample | L 102 -> 91, M 117 -> 113 |
+| `layer_m.comp`, `compile_m.py`: L's expansion stage at two windows per wave (`M_DEFS_DEC0L=-DWPW=1` for one) | expansion | 383 -> 366 |
+
+The network layers are separated by an execution-only barrier (`emptyStage`). A memory barrier there gives the same
+output and costs about 60 us (M 1520 -> 1592), so it stays.
+
 ## Preset L (runtime and shim path exist; experimental)
 
 Preset L (render preset 12, NVIDIA's model for Ultra Performance) is M's pipeline with its own input, expansion and

engine/compile_m.py

@@ -32,6 +32,8 @@ from test_m import pack_layer, defines
 EXACT = os.environ.get('M_EXACT') == '1'
 FAST = [] if EXACT or '-DM_RELAXED' in os.environ.get('M_DEFS', '').split() else ['-DM_RELAXED']
 FAST_WPW = [] if EXACT or os.environ.get('M_WPW') or 'WPW=' in os.environ.get('M_DEFS', '') else ['-DWPW=2']
+# preset L's expansion stage (dec0l) two windows per wave as well; its input stage (enc0l) is slower that way
+L_WPW = [] if EXACT or 'WPW=' in os.environ.get('M_DEFS', '') + os.environ.get('M_DEFS_DEC0L', '') else ['-DWPW=2']
 SHADERS = {'enc1': 'm_enc1', 'enc2': 'm_enc2', 'enc3_tube': 'm_tube', 'dec2': 'm_dec2', 'dec1': 'm_dec1'}
 
 
@@ -159,7 +161,7 @@ def compile_l(a, one, R):
     for name, out in SHADERS.items(): glsl('layer_m.comp', out + '.spv', defines(name) + ['-DTEXL=1'] + FAST + FAST_WPW)
     for planes, suffix in ((0, ''), (1, '_planes')):
         glsl('layer_m.comp', f'enc0_l{suffix}.spv', defines('enc0l') + ldr + [f'-DPLANES={planes}'] + FAST)
-        glsl('layer_m.comp', f'dec0_l{suffix}.spv', defines('dec0l') + [f'-DPLANES={planes}'] + FAST)
+        glsl('layer_m.comp', f'dec0_l{suffix}.spv', defines('dec0l') + [f'-DPLANES={planes}'] + FAST + L_WPW)
     # the reconstruction one render texel wider per tier than M's (4x4 taps): 11 (reconstruction/render >= 3),
     # 16 (the whole post_3_2 range), 12 (the compact tier: render/reconstruction <= 0.45, jitter within half a texel)
     for out, tile in (('post_l', 11), ('post_l_q', 16), ('post_l_c', 12)):

engine/compile_rr2.py

@@ -12,7 +12,7 @@ D4R_CUDA_REPLAY_DUMP_FILTER=rr2 D4R_CUDA_LAUNCH_STATS_FILTER=rr2 and preset hint
                        table of its PTX module; enc0-aux.bin is encoder 0's second table (argument 104)
   <stage>.args         the argument block of the capture's frame, pointers and handles cleared: the runtime keeps its
                        constants and fills what depends on the frame or the geometry (engine/runtime_rr2.inc)
-  exposure_f0.spv, exposure_f1.spv   engine/rr/exposure_f.comp
+  exposure_f0.spv, exposure_f1.spv   engine/rr/exposure_f.comp (exposure_f2.spv: both in one, 64-bit atomics)
   rr.txt               the manifest RModel reads ("D4RR0002")
 With --native, the stages engine/rr/rr2_m.comp implements (enc1 ... enc5, dec4 ... dec1, dec0's first kernel) are that
 shader with weights read from the DLL itself (engine/rr/rr2_pack.py) instead of the translation: bit-identical to
@@ -249,7 +249,7 @@ def main():
         if hybrid and split:   # the same block, bindings and regions; "rN": N words per window record, "s3": in three ranges
             lines.append(f'stage enc0g enc0g.spv 32 - {stage}.args {tex} {surf} {buf} {",".join(ptrs) or "-"} r{record} s3')
         print(f'{stage}: {kernel}, {len(args)} argument bytes' + (' (native block)' if hybrid else ''))
-    for p in (0, 1):
+    for p in (0, 1, 2):
         subprocess.run(['glslangValidator', '-V', '--target-env', 'vulkan1.3', f'-DPASS={p}',
                         str(HERE / 'rr' / 'exposure_f.comp'), '-o', str(out / f'exposure_f{p}.spv')],
                        check=True, stdout=subprocess.DEVNULL)

engine/down_m.comp

@@ -4,6 +4,7 @@
 // five-tap filter whose weights depend on the output pixel's parity, limits each pass to the range of its taps
 // (no ringing), and undoes the range compression and the exposure.
 #extension GL_EXT_shader_explicit_arithmetic_types : require
+#extension GL_EXT_control_flow_attributes : require
 #extension GL_GOOGLE_include_directive : require
 #include "output_store.glsl"
 // LDR (NGX's *_ldr_* kernel): the packed image is the display-referred colour itself.
@@ -13,8 +14,7 @@
 // D4RM0003: the runtime dispatches 8x64 output pixels per workgroup.
 #define BLOCK_X 8
 #define BLOCK_Y 64
-#define TILE_W (3 * BLOCK_X / 2 + 4)
-#define TILE_H (3 * BLOCK_Y / 2 + 4)
+#define TILE_H (3 * BLOCK_Y / 2 + 4)     // source rows of the block
 layout(local_size_x = BLOCK_X, local_size_y = BLOCK_Y) in;
 layout(set = 0, binding = 0, std430) readonly buffer Params
 {
@@ -27,8 +27,19 @@ layout(set = 0, binding = 0, std430) readonly buffer Params
 layout(set = 0, binding = 1) uniform sampler2D colorTex;        // render resolution; its alpha is passed through
 layout(set = 0, binding = 2, FINAL_OUTPUT_FORMAT) uniform writeonly image2D outImg;
 layout(set = 0, binding = 3, std430) readonly buffer Src { uint w[]; } src;
-shared f16vec4 tile[TILE_W * TILE_H];       // source texels around the block
 shared f16vec4 rows[TILE_H * BLOCK_X];      // after the horizontal pass
+// The horizontal pass reads the source texels straight into registers, without a shared tile: an invocation filters
+// SEG output columns of one source row from the 3 SEG / 2 + 4 texels around them.
+#define SEG 4
+#define SPAN (3 * SEG / 2 + 4)
+f16vec3 unpack(uint w)
+{
+#if LDR
+    return f16vec3(vec3(float(w & 2047u) * (1.0 / 2047.0), float((w >> 11) & 2047u) * (1.0 / 2047.0), float(w >> 22) * (1.0 / 1023.0)));
+#else
+    return f16vec3(vec3(float(w & 2047u) * 0.0011968734906986356, float((w >> 11) & 2047u) * 0.0011968734906986356, float(w >> 22) * 0.002394916955381632));
+#endif
+}
 // weights of taps -2 .. +2 for an even output coordinate; an odd one uses them mirrored
 // as binary16 bit patterns: NGX rounds its f32 constants to the nearest binary16 (cvt.rn), while glslang folds
 // float16_t(<literal>) toward zero (W0 became -0.041321 instead of -0.041351: 6% of the captured frame's outputs
@@ -49,25 +60,27 @@ void main()
 {
     const ivec2 lid = ivec2(gl_LocalInvocationID.xy), pix = ivec2(gl_GlobalInvocationID.xy);
     const ivec2 base = ivec2(3 * BLOCK_X / 2, 3 * BLOCK_Y / 2) * ivec2(gl_WorkGroupID.xy) - 2;
-    for (int i = lid.y * BLOCK_X + lid.x; i < TILE_W * TILE_H; i += BLOCK_X * BLOCK_Y)
+    // the horizontal pass, SEG output columns x0 .. x0 + SEG - 1 of source row y per item
+    const int zero = int(gl_NumWorkGroups.z) - 1;
+    for (int item = lid.y * BLOCK_X + lid.x; item < TILE_H * (BLOCK_X / SEG); item += BLOCK_X * BLOCK_Y)
     {
-        const ivec2 at = clamp(base + ivec2(i % TILE_W, i / TILE_W), ivec2(0), p.srcSize - 1);
-        const uint w = src.w[at.y * p.srcSize.x + at.x];
-#if LDR
-        tile[i] = f16vec4(vec4(float(w & 2047u) * (1.0 / 2047.0), float((w >> 11) & 2047u) * (1.0 / 2047.0), float(w >> 22) * (1.0 / 1023.0), 0.0));
-#else
-        tile[i] = f16vec4(vec4(float(w & 2047u) * 0.0011968734906986356, float((w >> 11) & 2047u) * 0.0011968734906986356, float(w >> 22) * 0.002394916955381632, 0.0));
-#endif
+        const int y = item / (BLOCK_X / SEG), x0 = item % (BLOCK_X / SEG) * SEG;
+        const int row = clamp(base.y + y, 0, p.srcSize.y - 1) * p.srcSize.x;
+        f16vec3 t[SPAN];
+        [[unroll]] for (int k = 0; k < SPAN; ++k)
+            t[k] = unpack(src.w[row + clamp(base.x + 3 * x0 / 2 + k, 0, p.srcSize.x - 1)]);
+        [[unroll]] for (int q = 0; q < SEG; ++q)
+        {
+            // column x0 + q (as even as q) is centred on texel 3 (q / 2) + 2 (q & 1) + 2 of the item's. Its parity stays a
+            // run-time value: with constant weights the compiler shares products between neighbouring filters and
+            // splits their fused multiply-adds (30% of the outputs change)
+            const int o = 3 * (q >> 1) + 2 * (q & 1) + 2;
+            rows[y * BLOCK_X + x0 + q].rgb = filter5(t[o - 2], t[o - 1], t[o], t[o + 1], t[o + 2], ((q + zero) & 1) != 0);
+        }
     }
     memoryBarrierShared(); barrier();
     const ivec2 c = ivec2(floor((vec2(pix) + 0.5) * 1.5)) - base;
-    for (int y = lid.y; y < TILE_H; y += BLOCK_Y)
-    {
-        const int o = y * TILE_W + c.x;
-        rows[y * BLOCK_X + lid.x].rgb = filter5(tile[o - 2].rgb, tile[o - 1].rgb, tile[o].rgb, tile[o + 1].rgb, tile[o + 2].rgb, (pix.x & 1) != 0);
-    }
-    memoryBarrierShared(); barrier();
-    // After both shared barriers: borrowed images may extend beyond the feature's true output rectangle.
+    // After the shared barrier: borrowed images may extend beyond the feature's true output rectangle.
     // M accepts even output dimensions and srcSize = outputSize/2*3, so this bounds check is exact.
     if (3 * pix.x >= 2 * p.srcSize.x || 3 * pix.y >= 2 * p.srcSize.y) return;
     const int o = c.y * BLOCK_X + lid.x;

engine/layer_m.comp

@@ -64,8 +64,8 @@
 #ifndef WPW
 #define WPW 1
 #endif
-#if WPW > 1 && KIND >= 3
-#error layer_m.comp: L's input and expansion stages stage their pixels per thread (WPW = 1)
+#if WPW > 2 && KIND >= 3
+#error layer_m.comp: L's input and expansion stages handle their 64 pixels one per invocation (WPW <= 2)
 #endif
 layout(local_size_x = 128 / WPW) in;
 // o: 0 g1, 1 bo, 2 g2, 3 b2 (f16); per head h at 4 + 4h: q, v (FP8), bias (f16), wo (FP8);
@@ -213,6 +213,15 @@ HAC bias16(uint off, uint stride) { HAC b; coopMatLoad(b, w16.w, off / 4u, strid
 // straight from f32 into FP8, two elements per saturating conversion; the activation in f32 as well
 #define SAT2_(dst, i, v) { fe4m3vec2 t_; saturatedConvertEXT(t_, (v)); (dst)[i] = t_.x; (dst)[(i) + 1] = t_.y; }
 #define MK8_F8(dst, a) { SAT2_(dst, 0, vec2((a)[0], (a)[1])) SAT2_(dst, 2, vec2((a)[2], (a)[3])) SAT2_(dst, 4, vec2((a)[4], (a)[5])) SAT2_(dst, 6, vec2((a)[6], (a)[7])) }
+// SATF2_: the same for a value that is never infinite, written as "an infinity becomes NaN, then saturate", which RADV
+// turns into the hardware conversion with FP16_OVFL instead of a NaN-preserving clamp per element and the conversion
+// (f2e4m3fn_satfn; the same result for every finite value and NaN, engine/rr/rr_exact.glsl FINITE). Used for the
+// attention's f32 sums of FP8 products (MK8_F8F) and the normalised tokens: the norm sums the squares of the very
+// values it scales, so an infinite channel makes the reciprocal 0 and the product NaN, never infinite. The binary16
+// accumulators (q, the MLP, the expansion, the residual) can overflow to infinity and keep SAT2_.
+vec2 finite2_(vec2 v) { return vec2(isinf(v.x) ? uintBitsToFloat(0x7fc00000u) : v.x, isinf(v.y) ? uintBitsToFloat(0x7fc00000u) : v.y); }
+#define SATF2_(dst, i, v) { fe4m3vec2 t_; saturatedConvertEXT(t_, finite2_(v)); (dst)[i] = t_.x; (dst)[(i) + 1] = t_.y; }
+#define MK8_F8F(dst, a) { SATF2_(dst, 0, vec2((a)[0], (a)[1])) SATF2_(dst, 2, vec2((a)[2], (a)[3])) SATF2_(dst, 4, vec2((a)[4], (a)[5])) SATF2_(dst, 6, vec2((a)[6], (a)[7])) }
 #define GELU2_(dst, i, h) { const vec2 x_ = vec2((h)[i], (h)[(i) + 1]); const vec2 y_ = clamp(x_, vec2(-2.0), vec2(2.0)); \
     const vec2 g_ = fma(y_, fma(-abs(y_), vec2(0.08108133), vec2(0.41216981)), vec2(0.5)); SAT2_(dst, i, x_ * g_) }
 #define GELU8(dst, h) { GELU2_(dst, 0, h) GELU2_(dst, 2, h) GELU2_(dst, 4, h) GELU2_(dst, 6, h) }
@@ -222,12 +231,13 @@ HAC bias16(uint off, uint stride) { HAC b; coopMatLoad(b, w16.w, off / 4u, strid
     const float16_t q4_ = CL8(Q16((a)[4])), q5_ = CL8(Q16((a)[5])), q6_ = CL8(Q16((a)[6])), q7_ = CL8(Q16((a)[7])); \
     (dst)[0] = floate4m3_t(q0_); (dst)[1] = floate4m3_t(q1_); (dst)[2] = floate4m3_t(q2_); (dst)[3] = floate4m3_t(q3_); \
     (dst)[4] = floate4m3_t(q4_); (dst)[5] = floate4m3_t(q5_); (dst)[6] = floate4m3_t(q6_); (dst)[7] = floate4m3_t(q7_); }
+#define MK8_F8F MK8_F8
 #endif
 // the normalised tokens into FP8 (inside NM1 / N2: x0..x7 the token, g0 / g1 the gains)
 #ifdef M_RELAXED
 #define H1N(N, r_) { const float f_ = float(r_); \
-        SAT2_(h1[N], 0, vec2(x0, x1) * (f_ * vec2(g0.xy))) SAT2_(h1[N], 2, vec2(x2, x3) * (f_ * vec2(g0.zw))) \
-        SAT2_(h1[N], 4, vec2(x4, x5) * (f_ * vec2(g1.xy))) SAT2_(h1[N], 6, vec2(x6, x7) * (f_ * vec2(g1.zw))) }
+        SATF2_(h1[N], 0, vec2(x0, x1) * (f_ * vec2(g0.xy))) SATF2_(h1[N], 2, vec2(x2, x3) * (f_ * vec2(g0.zw))) \
+        SATF2_(h1[N], 4, vec2(x4, x5) * (f_ * vec2(g1.xy))) SATF2_(h1[N], 6, vec2(x6, x7) * (f_ * vec2(g1.zw))) }
 #else
 #define H1N(N, r_) { \
         h1[N][0] = floate4m3_t(CL8(x0 * float16_t(r_ * g0.x))); h1[N][1] = floate4m3_t(CL8(x1 * float16_t(r_ * g0.y))); \
@@ -356,7 +366,6 @@ float16_t gelu1(float16_t x)
     precise f16vec2 C_ = f16vec2(EV(T, 4), EV(T, 5)), D_ = f16vec2(EV(T, 6), EV(T, 7)); \
     P0_[T] = (A_ * A_) + (B_ * B_); P1_[T] = (C_ * C_) + (D_ * D_); }
 
-shared f16vec4 VSH[4 * 2 * 16 * VWD / 4];   // per wave and head: the V tiles, transposed through shared memory
 #define MERGE (KIND == 1 || KIND == 3)
 #define EXPAND (KIND == 2 || KIND == 4)
 #if MERGE
@@ -375,6 +384,34 @@ shared float16_t X0S[64 * XS];              // the embedding's output x0, the bl
 const uint HS = 40u;
 shared float16_t HO[64 * HS];               // the block's f16 output, window-major, for the per-pixel head
 #endif
+// VSH: per wave and head the V tiles, transposed through shared memory (4 KiB, 512 bytes a tile). It is only accessed
+// as matrices, and only during the attention, so it takes the place of a buffer that is dead then: the merge rows
+// (written after the attention), the expansion (read before it) or, for KIND 4, the head's input (written after it).
+// Shared memory caps how many waves this GPU runs, and ACO gives a shader the registers that number of waves leaves.
+#if MERGE
+#define VSH NB
+#define VSH_UNIT 1u
+#if 32 * (4 * C + 8) < 4096
+#error layer_m.comp: the V tiles no longer fit in NB
+#endif
+#elif KIND == 2
+#define VSH NEX
+#define VSH_UNIT 1u
+#if 16 * (4 * C + 8) < 4096
+#error layer_m.comp: the V tiles no longer fit in NEX
+#endif
+#elif KIND == 4
+#define VSH HO
+#define VSH_UNIT 2u
+#if 64 * 40 * 2 < 4096
+#error layer_m.comp: the V tiles no longer fit in HO
+#endif
+#else
+shared floate4m3_t VSH[4 * 2 * 16 * VWD * 2];
+#define VSH_UNIT 1u
+#endif
+#define VSH_AT(t) ((t) * (512u / VSH_UNIT))      // tile t, in the buffer's elements
+#define VSH_STRIDE (VWD * 2u / VSH_UNIT)          // a row of VWD binary16
 #define SYNC memoryBarrierShared(); barrier();
 #if KIND == 4
 // NVIDIA's f32 activations of the head (rrlite_dec0_4x4's output code): a sigmoid rounded to binary16, and a
@@ -459,8 +496,7 @@ void main()
 #endif
 #if KIND == 3
     // ------------------------------------------------------------ L's input: features, embedding -> x0
-    FORW {
-        WINV
+    {
         const float e1 = ep.expo == 0.0 ? 1.0 : ep.expo, e2 = ep.expoPrev == 0.0 ? 1.0 : ep.expoPrev;
         const uint li = gl_LocalInvocationID.x;
         if (li < 64u)
@@ -473,7 +509,10 @@ void main()
             const uint at = (((py >> 2) * 2u + (px >> 2)) * 16u + (py & 3u) * 4u + (px & 3u)) * FS;
             for (uint k = 0u; k < 32u; ++k) FEW[at + k] = fe[k];
         }
+    }
     SYNC
+    FORW {
+        WINV
         // the embedding has the shape of an MLP chunk: 32 features -> 32 hidden (ReLU, FP8) -> C
         B8 f0, f1;
         coopMatLoad(f0, FEW, mt * 16u * FS, FS, CM);
@@ -495,8 +534,8 @@ void main()
         for (int i = 0; i < 8; ++i) X0S[xb + 16u * (N) + uint(i)] = float16_t(x_[i]); }
         EACH_T(EMB_N)
 #undef EMB_N
-        SYNC
     }
+    SYNC
 #endif
     // ------------------------------------------------------------ norm 1: sum of squares, residual and h1
     FORW {
@@ -542,6 +581,9 @@ void main()
     H1N(N, rn) }
     EACH_T(NM1)
     }
+#if KIND == 2
+    SYNC   // every wave has read the block's input from NEX: the V tiles take its place
+#endif
     // ------------------------------------------------------------ attention, one head at a time
     for (uint hd = 0u; hd < uint(NH); ++hd)
     {
@@ -593,8 +635,8 @@ void main()
         // V of every window goes to its shared buffer now, so it does not occupy registers during the scores
         FORW {
             const uint mt = mt0 + w;
-            coopMatStore(v[w][0], VSH, (mt * 2u + 0u) * 64u, VWD / 4u, RM);
-            coopMatStore(v[w][1], VSH, (mt * 2u + 1u) * 64u, VWD / 4u, RM);
+            coopMatStore(v[w][0], VSH, VSH_AT(mt * 2u + 0u), VSH_STRIDE, RM);
+            coopMatStore(v[w][1], VSH, VSH_AT(mt * 2u + 1u), VSH_STRIDE, RM);
         }
         }
         // scores: q8(Q) against itself, rows are the queries and columns the keys, plus the relative position bias
@@ -624,8 +666,8 @@ void main()
         // V through this wave's shared buffer, then O^T = V^T x P^T: rows are the head channels, columns the tokens
         ACC o[2];
         A16 vt0, vt1;
-        coopMatLoad(vt0, VSH, (mt * 2u + 0u) * 64u, VWD / 4u, RM);
-        coopMatLoad(vt1, VSH, (mt * 2u + 1u) * 64u, VWD / 4u, RM);
+        coopMatLoad(vt0, VSH, VSH_AT(mt * 2u + 0u), VSH_STRIDE, RM);
+        coopMatLoad(vt1, VSH, VSH_AT(mt * 2u + 1u), VSH_STRIDE, RM);
 #if TEXL
         // P V is the block's one product of two f16 operands, and its f32 sum over the 16 keys is not exact: ZLUDA's
         // lowering of m16n8k16 puts key pair (j & 1) * 4 + 2 h + (j >> 1) in K dword j of half h (swin_tok.h
@@ -635,7 +677,7 @@ void main()
         o[0] = coopMatMulAdd(vt0, pb, ACC(0.0));
         o[1] = coopMatMulAdd(vt1, pb, ACC(0.0));
         R16(o[0]); R16(o[1]);
-        MK8_F8(OB_[w][0], o[0]); MK8_F8(OB_[w][1], o[1]);
+        MK8_F8F(OB_[w][0], o[0]); MK8_F8F(OB_[w][1], o[1]);
         }
         // the output projection of this head's 32 attention channels, straight into the residual stream
 #define WO_N(N) { \
@@ -722,12 +764,15 @@ void main()
     }
 #if KIND == 4
     // ------------------------------------------------------------ L's expansion head: per-pixel reconstruction inputs
+    SYNC   // every wave is done with the V tiles in HO
     FORW {
         WINV
         const uint hb = (mt * 16u + row) * HS + ch8;
         for (uint nt = 0u; nt < NT; ++nt)
             for (int i = 0; i < 8; ++i) HO[hb + 16u * nt + uint(i)] = float16_t(res[nt][i]);
+    }
     SYNC
+    {
         const uint li = gl_LocalInvocationID.x;
         const uint px = li & 7u, py = li >> 3;
         const ivec2 at = ivec2(X0 + int(px), Y0 + int(py));
@@ -763,6 +808,7 @@ void main()
 #endif
     // ------------------------------------------------------------ patch merge: 2x2 tokens (4C) -> XC, half resolution
 #if MERGE
+    SYNC   // every wave is done with the V tiles in NB
     FORW {
         WINV
         // the 16 merged rows are gathered from the four waves: the row is the merged token of the lane's pixel

engine/post_m.comp

@@ -154,6 +154,42 @@ void main()
     // three render texel rows and its own-cell residual rounding exactly as in the 16x16 shader
     const int y0 = 2 * gid.y;
     const bool inside0 = gid.x < p.hiSize.x && y0 < p.hiSize.y;
+    const vec2 pos0 = fma(vec2(gid.x, y0) + 0.5, p.hi2rnd, p.rndOffs), pos1 = fma(vec2(gid.x, y0 + 1) + 0.5, p.hi2rnd, p.rndOffs);
+    // The Gaussians of both rows at once: they have the same column (dx) and texels and differ only in dy, so one set
+    // of tile reads serves both, and each row's operations are the ones it had on its own, lane by lane (dy, the
+    // exponents and the third accumulator are binary16 pairs, element 0 the even row).
+    f16vec4 a0[2] = {f16vec4(0.0), f16vec4(0.0)}, a1[2] = {f16vec4(0.0), f16vec4(0.0)};
+    f16vec2 a2 = f16vec2(0.0);
+    float W[2] = {0.0, 0.0};
+    if (inside0)
+    {
+        // NVIDIA's kernel handles rows in pairs and gives the odd row the even row's three texel rows
+        const ivec2 base = ivec2(floor(pos0 - 0.5));
+        const f16vec2 d0 = f16vec2((vec2(base) + 0.5) - pos0);
+        // Native paired rows round the odd row's own-cell residual first, then adjust it in f16.
+        // Casting the even-cell residual once changes Gaussian weights at half-rounding boundaries.
+        const int ownBaseY = int(floor(pos1.y - 0.5));
+        const float16_t d1y = ownBaseY != base.y ? float16_t((float(ownBaseY) + 0.5) - pos1.y) - float16_t(1.0)
+                                                 : float16_t((float(base.y) + 0.5) - pos1.y);
+        const f16vec2 dy0 = f16vec2(d0.y, d1y);
+        for (int i = TAP_LO; i <= TAP_HI; ++i)
+        {
+            const float16_t dx = d0.x + float16_t(i), sq = dx * dx;
+            for (int j = TAP_LO; j <= TAP_HI; ++j)
+            {
+                const int k = (base.y - origin.y + j) * TILE_W + (base.x - origin.x + i);
+                const f16vec4 t0 = tile0[k], t1 = tile1[k], t2 = tile2[k];
+                const f16vec2 dy = dy0 + float16_t(j);
+                precise f16vec2 e = (f16vec2(sq * t2.y) + (dy * dy) * t2.z) + (dx * dy) * t2.w;
+                const float w0 = exp2(float(e.x)), w1 = exp2(float(e.y));
+                const float16_t wh0 = float16_t(w0), wh1 = float16_t(w1);
+                W[0] += w0; W[1] += w1;
+                a0[0] = fma(f16vec4(wh0), t0, a0[0]); a1[0] = fma(f16vec4(wh0), t1, a1[0]);
+                a0[1] = fma(f16vec4(wh1), t0, a0[1]); a1[1] = fma(f16vec4(wh1), t1, a1[1]);
+                a2 = fma(f16vec2(wh0, wh1), f16vec2(t2.x), a2);
+            }
+        }
+    }
     vec3 outA0 = vec3(0.0), outA1 = vec3(0.0);
     for (int r = 0; r < 2; ++r)
     {
@@ -163,39 +199,10 @@ void main()
         if (inside)
         {
             const vec2 centre = vec2(pix) + 0.5;
-            const vec2 pos = fma(centre, p.hi2rnd, p.rndOffs);
-            // NVIDIA's kernel handles rows in pairs and gives the odd row the even row's three texel rows
-            const float evenY = fma(float(pix.y & ~1) + 0.5, p.hi2rnd.y, p.rndOffs.y);
-            const ivec2 near = ivec2(floor(pos)), base = ivec2(floor(pos.x - 0.5), floor(evenY - 0.5));
-            f16vec2 d0 = f16vec2((vec2(base) + 0.5) - pos);
-            if ((pix.y & 1) != 0)
-            {
-                // Native paired rows round the odd row's own-cell residual first, then adjust it in f16.
-                // Casting the even-cell residual once changes Gaussian weights at half-rounding boundaries.
-                const int ownBaseY = int(floor(pos.y - 0.5));
-                if (ownBaseY != base.y)
-                    d0.y = float16_t((float(ownBaseY) + 0.5) - pos.y) - float16_t(1.0);
-            }
-            f16vec4 a0 = f16vec4(0.0), a1 = f16vec4(0.0);
-            float16_t a2 = float16_t(0.0);
-            float W = 0.0;
-            for (int i = TAP_LO; i <= TAP_HI; ++i)
-            {
-                const float16_t dx = d0.x + float16_t(i), sq = dx * dx;
-                for (int j = TAP_LO; j <= TAP_HI; ++j)
-                {
-                    const int k = (base.y - origin.y + j) * TILE_W + (base.x - origin.x + i);
-                    const f16vec4 t2 = tile2[k];
-                    const float16_t dy = d0.y + float16_t(j);
-                    precise float16_t e = (sq * t2.y + (dy * dy) * t2.z) + (dx * dy) * t2.w;
-                    const float w = exp2(float(e));
-                    const float16_t wh = float16_t(w);
-                    W += w;
-                    a0 = fma(f16vec4(wh), tile0[k], a0); a1 = fma(f16vec4(wh), tile1[k], a1); a2 = fma(wh, t2.x, a2);
-                }
-            }
-            const float inv = 1.0 / W;
-            const vec3 cur = inv * vec3(a0.rgb), m = inv * vec3(a0.w, a1.x, a1.y), g = inv * vec3(a1.z, a1.w, a2);
+            const vec2 pos = r == 0 ? pos0 : pos1;
+            const ivec2 near = ivec2(floor(pos));
+            const float inv = 1.0 / W[r];
+            const vec3 cur = inv * vec3(a0[r].rgb), m = inv * vec3(a0[r].w, a1[r].x, a1[r].y), g = inv * vec3(a1[r].z, a1[r].w, a2[r]);
             // the two most likely motion candidates
             const vec2 mvPos = any(notEqual(textureSize(mvCTex, 0), p.rndSize))
                 ? clamp(pos, vec2(0.5), vec2(p.rndSize) - 0.5) : pos;
@@ -214,6 +221,9 @@ void main()
             // the fallback colour is only fetched where the history sample leaves the image or the history was reset:
             // the value is the same either way, so computing it in the branch saves the fetch on ordinary pixels
             const bool out1 = q1.x < 0.0 || q1.y < 0.0 || q1.x > size.x || q1.y > size.y || p.reset != 0u;
+            // the second candidate often has the first one's vector (where the motion is uniform): the same position,
+            // so the first one's sample or fallback (a NaN is never equal, and +0 and -0 give the same position)
+            const bool same = all(equal(second.xy, big.xy));
 #if PRESET_L
             // linear light: a sample's expansion (when the previous exposure is positive) or the colour texel
             const vec3 h1 = out1 ? texelFetch(colorTex, colourTexel(near), 0).rgb : historyLinear(q1, p.expoPrev == 0.0 ? 1.0 : p.expoPrev);
@@ -221,7 +231,9 @@ void main()
             if (b2 > 0.0)
             {
                 const vec2 q2 = centre + second.xy;
-                if (q2.x < 0.0 || q2.y < 0.0 || q2.x > size.x || q2.y > size.y || p.reset != 0u)
+                if (same)
+                    h2 = h1;
+                else if (q2.x < 0.0 || q2.y < 0.0 || q2.x > size.x || q2.y > size.y || p.reset != 0u)
                     h2 = out1 ? h1 : texelFetch(colorTex, colourTexel(near), 0).rgb;
                 else h2 = historyLinear(q2, p.expoPrev == 0.0 ? 1.0 : p.expoPrev);
             }
@@ -232,14 +244,16 @@ void main()
             if (b2 > 0.0)
             {
                 const vec2 q2 = centre + second.xy;
-                if (q2.x < 0.0 || q2.y < 0.0 || q2.x > size.x || q2.y > size.y || p.reset != 0u)
+                if (same)
+                    h2 = h1;
+                else if (q2.x < 0.0 || q2.y < 0.0 || q2.x > size.x || q2.y > size.y || p.reset != 0u)
                     h2 = out1 ? h1 : compress(e2 * texelFetch(colorTex, colourTexel(near), 0).rgb);
                 else h2 = history(q2);
             }
             const vec3 hist = e1 * (expand(fma(vec3(b1), h1, b2 * h2)) * (1.0 / e2));
 #endif
             const vec3 lin = expand(cur);
-            const float g0 = clamp(p.gain * W * g.x, 0.0, 1.0);
+            const float g0 = clamp(p.gain * W[r] * g.x, 0.0, 1.0);
             const float ga = fma(g0, 1.0 - g.y, g.y), gb = fma(g0, 1.0 - g.z, g.z);
             outA = fma(lin - hist, vec3(ga), hist);
             if (isinf(outA.x + outA.y + outA.z) || isnan(outA.x + outA.y + outA.z)) outA = lin;

engine/replay_m.cpp

@@ -107,7 +107,7 @@ int main(int argc, char** argv) try {
     CK(vkCreateCommandPool(c.device, &pi, nullptr, &c.pool));
     VkCommandBufferAllocateInfo cai{VK_STRUCTURE_TYPE_COMMAND_BUFFER_ALLOCATE_INFO}; cai.commandPool = c.pool; cai.level = VK_COMMAND_BUFFER_LEVEL_PRIMARY; cai.commandBufferCount = 1;
     VkCommandBuffer cb; CK(vkAllocateCommandBuffers(c.device, &cai, &cb));
-    VkQueryPoolCreateInfo qpi{VK_STRUCTURE_TYPE_QUERY_POOL_CREATE_INFO}; qpi.queryType = VK_QUERY_TYPE_TIMESTAMP; qpi.queryCount = 7;
+    VkQueryPoolCreateInfo qpi{VK_STRUCTURE_TYPE_QUERY_POOL_CREATE_INFO}; qpi.queryType = VK_QUERY_TYPE_TIMESTAMP; qpi.queryCount = 17;
     VkQueryPool qp; CK(vkCreateQueryPool(c.device, &qpi, nullptr, &qp));
     VkPhysicalDeviceProperties props; vkGetPhysicalDeviceProperties(c.physical, &props);
     auto submit = [&] {
@@ -151,7 +151,7 @@ int main(int argc, char** argv) try {
                 if (data.size() != sizes[k]) throw std::runtime_error("wrong input size");
                 frameInputs.push_back(std::move(data));
             }
-        std::vector<double> times, walls; std::vector<std::vector<double>> stage(6);
+        std::vector<double> times, walls; std::vector<std::vector<double>> stage(6), layer(10);
         const auto now = [] { return std::chrono::duration<double, std::micro>(std::chrono::steady_clock::now().time_since_epoch()).count(); };
         // D4R_M_RESET_ALL=1: every frame is a reset frame, as the D3D12 harness passes them without its temporal option
         const bool resetAll = std::getenv("D4R_M_RESET_ALL") != nullptr;
@@ -169,7 +169,7 @@ int main(int argc, char** argv) try {
                 s.renderWidth = rw; s.renderHeight = rh;
                 const double t0 = now();
                 engine.setFrame({color.view, motion.view, depth.view}, s);
-                vkCmdResetQueryPool(cb, qp, 0, 7);
+                vkCmdResetQueryPool(cb, qp, 0, 17);
                 double networkUs = 0;
                 if (hip) {
 #ifdef D4R_ENGINE_HIP
@@ -186,12 +186,16 @@ int main(int argc, char** argv) try {
                 if (pass == 0) vkCmdCopyImageToBuffer(cb, engine.outputImage(), VK_IMAGE_LAYOUT_GENERAL, out.buffer, 1, &region);
                 submit();
                 if (pass) walls.push_back(now() - t0);
-                uint64_t t[7]; CK(vkGetQueryPoolResults(c.device, qp, 0, 7, sizeof(t), t, 8, VK_QUERY_RESULT_64_BIT | VK_QUERY_RESULT_WAIT_BIT));
+                uint64_t t[17] = {};
+                // D4R_ENGINE_M_LAYER_STAMPS=1: the ten network layers' stamps as well (queries 7 - 16)
+                static const bool layerStamps = std::getenv("D4R_ENGINE_M_LAYER_STAMPS") != nullptr && !hip;
+                CK(vkGetQueryPoolResults(c.device, qp, 0, layerStamps ? 17 : 7, sizeof(t), t, 8, VK_QUERY_RESULT_64_BIT | VK_QUERY_RESULT_WAIT_BIT));
                 const double tick = props.limits.timestampPeriod / 1000.0;
                 double d[6]; for (int k = 0; k < 6; ++k) d[k] = double(t[k + 1] - t[k]) * tick;
                 if (hip) d[2] = networkUs; // measured HIP network time; the complete span includes host handoffs
                 const double us = double(t[6] - t[0]) * tick;
                 if (pass) for (int k = 0; k < 6; ++k) stage[k].push_back(d[k]);
+                if (pass && layerStamps) for (int l = 0; l < 10; ++l) layer[l].push_back(double(t[7 + l] - (l ? t[6 + l] : t[2])) * tick);
                 if (pass == 0 && std::getenv("D4R_M_DUMP")) {
                     // development: every intermediate of this frame, to PREFIX-NAME-FRAME.bin
                     for (const auto& d : engine.debugResources()) {
@@ -216,6 +220,11 @@ int main(int argc, char** argv) try {
             std::sort(walls.begin(), walls.end()); std::printf("  wall clock, record to completion: median %.1f us, min %.1f us\n", walls[walls.size() / 2], walls[0]);
             const char* names[6] = {"exposure", "input", "network", "expansion", "reconstruction", "downsample"};
             for (int k = 0; k < 6; ++k) { std::sort(stage[k].begin(), stage[k].end()); std::printf("  %-14s median %.1f us\n", names[k], stage[k][stage[k].size() / 2]); }
+            if (!layer[0].empty()) {
+                std::printf("  layers (us):");
+                for (auto& v : layer) { std::sort(v.begin(), v.end()); std::printf(" %.1f", v[v.size() / 2]); }
+                std::printf("\n");
+            }
         }
         CK(vkDeviceWaitIdle(c.device));
     }

engine/replay_rr.cpp

@@ -158,6 +158,8 @@ int main(int argc, char** argv) try {
     if (std::getenv("D4R_RR_TEST_NO_ROBUST")) { f2.features.robustBufferAccess = VK_FALSE; v13.robustImageAccess = VK_FALSE; }
     // D4R_RR_TEST_NO_COOPMAT=1: a device like vkd3d-proton's, without the cooperative-matrix and float8 features
     if (std::getenv("D4R_RR_TEST_NO_COOPMAT")) { di.enabledExtensionCount = 0; v13.pNext = nullptr; }
+    // D4R_RR_TEST_NO_INT64_ATOMICS=1: no 64-bit buffer atomics (preset F's exposure then takes its two passes)
+    if (std::getenv("D4R_RR_TEST_NO_INT64_ATOMICS")) v12.shaderBufferInt64Atomics = v12.shaderSharedInt64Atomics = VK_FALSE;
     CK(vkCreateDevice(c.physical, &di, nullptr, &c.device));
     vkGetDeviceQueue(c.device, family, 0, &c.queue);
     VkCommandPoolCreateInfo pi{VK_STRUCTURE_TYPE_COMMAND_POOL_CREATE_INFO}; pi.queueFamilyIndex = family;
@@ -169,7 +171,7 @@ int main(int argc, char** argv) try {
     std::string key;
     uint32_t rw, rh, ow, oh;
     frames >> key >> rw >> rh >> ow >> oh;
-    d4r::RModel model(c.physical, c.device, modelDir);
+    d4r::RModel model(c.physical, c.device, modelDir, VK_NULL_HANDLE, v12.shaderBufferInt64Atomics == VK_TRUE);
     d4r::REngine engine(model, rw, rh, ow, oh);
     const char* names[] = {"color", "motion", "depth", "normals", "roughness", "diffuse", "specular"};
     const VkFormat formats[] = {VK_FORMAT_R16G16B16A16_SFLOAT, VK_FORMAT_R16G16_SFLOAT, VK_FORMAT_R32_SFLOAT,

engine/rr/README.md

@@ -615,3 +615,17 @@ was measured and not built for D / E (`build/rr/f/nv/dstage.sh`, `dnl.sh`, level
   loads between a window's waves) decoder 0 takes the same time (803-818 us), encoder 1 the same, decoder 1 10% less.
   These windows also have 96 key tokens (3 KiB each for the tokens, the keys and the values at C = 32), which does
   not fit the 4 KiB a single wave can hold at full rate.
+
+### Byte-identical speed-ups (2026-10-19)
+
+RX 9070 XT, RADV 26.2.2, 1280x720 -> 2560x1440 with `D4R_ENGINE_RR_STAGE_STAMPS=1`. Every image of every frame is
+byte-identical to the previous model on the replay and stress inputs, at levels 0 to 3.
+
+- `ptx2rr.py`: OCML's tanh (decoder 0's second kernel) in f32 where f64 was not needed. Exp's error term becomes
+  Dekker's exact product error; the polynomial below 0.625 stays in f64, because the driver folds away the f32
+  rounding between its f64 steps. Identical to the previous function on all 2^32 inputs. dec0_1 475 -> 446 us, frame
+  3.568 -> 3.537 ms.
+- `exposure_f.comp` PASS 2: the exposure in one pass as NGX does it, each block adding its fixed-point sum to a 64-bit
+  total with an atomic add and the last block finishing. It needs `shaderBufferInt64Atomics`, which vkd3d-proton
+  enables; the shim reads the device's features through `GetDeviceFeatures`. Without it, or with a model that has no
+  `exposure_f2.spv`, the two passes run (`D4R_RR_TEST_NO_INT64_ATOMICS=1` forces them). Exposure 25.2 -> 14.9 us.

engine/rr/exposure_f.comp

@@ -1,33 +1,39 @@
 #version 460
-// Ray Reconstruction preset F (rr2_*): the auto exposure of rr2_histogram_auto_exposure_basic_kernel, in two passes.
+// Ray Reconstruction preset F (rr2_*): the auto exposure of rr2_histogram_auto_exposure_basic_kernel.
 //
 // NVIDIA's kernel sums a clamped luma per 16x16 block (a pairwise tree over the block's threads), converts each block
 // sum to fixed point (x 2^32 / 10000, round to nearest) and adds them with a 64-bit atomic; the last block to finish
-// turns the total into the exposure. The integer total does not depend on the order of the blocks, so PASS 0 writes
-// every block's fixed-point sum and PASS 1 (one workgroup) adds them and does what the last block does:
+// turns the total into the exposure:
 //   exposure = 0.23547088 / max(mean luma, 0.005)               -> the 1x1 exposure image of this frame
 //   values   = { e, p, 1 / e, clamp(log2(e / p), -6, 6) / 6 }   -> the 16-byte buffer encoder 0 and the output passes read
 // with e the exposure truncated to binary16 and p the previous frame's (1 where either is zero).
+// PASS 2 is that kernel (64-bit atomics on a storage buffer: shaderBufferInt64Atomics). Without them the runtime runs
+// two passes: PASS 0 writes every block's fixed-point sum and PASS 1 (one workgroup) adds them and does what the last
+// block does. The total is an integer: it does not depend on the order of the blocks, so both give the same values.
 #extension GL_EXT_shader_explicit_arithmetic_types : require
 #extension GL_EXT_control_flow_attributes : require
 #ifndef PASS
 #define PASS 0
 #endif
+#if PASS == 2
+#extension GL_EXT_shader_atomic_int64 : require
+#extension GL_KHR_memory_scope_semantics : require
+#endif
 layout(set = 0, binding = 0, std430) readonly buffer Params {
     ivec2 last;        // render size - 1
     uint blocksX, blocks;
     float invCount;    // 10000 / (2^32 * render pixels)
     uint reset;        // the kernel's argument 32: no log ratio when set
 } p;
-#if PASS == 0
+#if PASS != 1
 // One thread per column of the 16x16 block. NVIDIA's tree over the block's 256 threads (index 16 y + x: i += i + 128,
 // then 64, 32, ... 1) is, per column, ((r0 + r8) + (r4 + r12)) + ((r2 + r10) + (r6 + r14)) of the even rows plus the
 // same of the odd rows, and then the columns (c + c8, + 4, + 2, + 1).
 layout(local_size_x = 16) in;
 layout(set = 0, binding = 1) uniform sampler2D color;
-layout(set = 0, binding = 2, std430) writeonly buffer Sums { uint64_t s[]; } sums;
 shared float col[16];
-void main()
+// the block's fixed-point sum, in thread 0
+uint64_t blockSum()
 {
     const uint tx = gl_LocalInvocationID.x;
     const int x = int(gl_GlobalInvocationID.x), y0 = int(gl_WorkGroupID.y) * 16;
@@ -49,33 +55,37 @@ void main()
     precise float c = w1 + w0;
     col[tx] = c;
     barrier();
-    if (tx == 0u) {
+    if (tx != 0u) return 0ul;
     precise float d[8];
     [[unroll]] for (int k = 0; k < 8; ++k) d[k] = col[k + 8] + col[k];
     precise float e0 = d[4] + d[0], e1 = d[5] + d[1], e2 = d[6] + d[2], e3 = d[7] + d[3];
     precise float f0 = e2 + e0, f1 = e3 + e1;
     precise float g = f1 + f0;
-        sums.s[gl_WorkGroupID.y * p.blocksX + gl_WorkGroupID.x] = uint64_t(roundEven(g * 429496.71875));
+    return uint64_t(roundEven(g * 429496.71875));
 }
+#endif
+#if PASS == 0
+layout(set = 0, binding = 2, std430) writeonly buffer Sums { uint64_t s[]; } sums;
+void main()
+{
+    const uint64_t s = blockSum();
+    if (gl_LocalInvocationID.x == 0u) sums.s[gl_WorkGroupID.y * p.blocksX + gl_WorkGroupID.x] = s;
 }
 #else
+#if PASS == 1
 layout(local_size_x = 64) in;
 shared uint64_t part[64];
 layout(set = 0, binding = 1, std430) readonly buffer Sums { uint64_t s[]; } sums;
-layout(set = 0, binding = 2, std430) writeonly buffer Values { vec4 v; } values;
-layout(set = 0, binding = 3) uniform sampler2D previous;
-layout(set = 0, binding = 4, r16f) uniform writeonly image2D current;
-void main()
+#else
+// s[0]: the total, s[1]: the blocks done; the last block leaves both 0 for the next frame
+layout(set = 0, binding = 2, std430) buffer Sums { uint64_t s[]; } sums;
+#endif
+layout(set = 0, binding = PASS + 1, std430) writeonly buffer Values { vec4 v; } values;
+layout(set = 0, binding = PASS + 2) uniform sampler2D previous;
+layout(set = 0, binding = PASS + 3, r16f) uniform writeonly image2D current;
+// the last block's work, from the frame's total
+void finish(uint64_t total)
 {
-    // the block sums are integers: any order. 64 strided partial sums, then thread 0 finishes
-    const uint tid = gl_LocalInvocationID.x;
-    uint64_t mine = 0ul;
-    for (uint b = tid; b < p.blocks; b += 64u) mine += sums.s[b];
-    part[tid] = mine;
-    barrier();
-    if (tid != 0u) return;
-    uint64_t total = 0ul;
-    for (uint b = 0u; b < 64u; ++b) total += part[b];
     precise float mean = p.invCount * float(total);
     precise float e = 0.23547088 / max(mean, 0.005);
     imageStore(current, ivec2(0), vec4(e));
@@ -90,4 +100,33 @@ void main()
     if (p.reset == 0u) ratio = clamp(log2(max(cur, 1e-8) / max(prev, 1e-8)), -6.0, 6.0) * 0.16666667;
     values.v = vec4(cur, prev, 1.0 / cur, ratio);
 }
+#if PASS == 1
+void main()
+{
+    // the block sums are integers: any order. 64 strided partial sums, then thread 0 finishes
+    const uint tid = gl_LocalInvocationID.x;
+    uint64_t mine = 0ul;
+    for (uint b = tid; b < p.blocks; b += 64u) mine += sums.s[b];
+    part[tid] = mine;
+    barrier();
+    if (tid != 0u) return;
+    uint64_t total = 0ul;
+    for (uint b = 0u; b < 64u; ++b) total += part[b];
+    finish(total);
+}
+#else
+void main()
+{
+    const uint64_t s = blockSum();
+    if (gl_LocalInvocationID.x != 0u) return;
+    atomicAdd(sums.s[0], s, gl_ScopeDevice, 0, gl_SemanticsRelaxed);
+    // the release orders this block's add before its count, and the last block's acquire every block's add before
+    // its read of the total
+    if (atomicAdd(sums.s[1], 1ul, gl_ScopeDevice, gl_StorageSemanticsBuffer, gl_SemanticsAcquireRelease) != uint64_t(p.blocks - 1u))
+        return;
+    const uint64_t total = atomicExchange(sums.s[0], 0ul, gl_ScopeDevice, 0, gl_SemanticsRelaxed);
+    atomicExchange(sums.s[1], 0ul, gl_ScopeDevice, 0, gl_SemanticsRelaxed);
+    finish(total);
+}
+#endif
 #endif

engine/rr/ptx2rr.py

@@ -121,31 +121,49 @@ void sust_p1(writeonly image2D s, int x, int y, float v) { sust_p4(s, x, y, vec4
 // UNORM 10:10:10:2 surfaces (the pooled histories) are RGBA32F images holding the quantised values, like ZLUDA's.
 // tanh.approx.f32 as the reference computes it: ROCm's OCML tanh (device-libs ocml/src/tanhF.cl) with LLVM's AMDGPU
 // lowering of exp (lowerFEXP, fast-FMA path): the same fused operations, hardware exp2, ldexp and reciprocal.
-// (a fused f32 a * b + c through f64: the product is exact there and the sum is rounded to 53 bits before the f32
-// rounding. The driver's own f32 fma is a product and a sum whenever its product also exists on its own, as in
-// exp's ph / fma(x, c, -ph) pair.)
-float fmad(float a, float b, float c) { precise double s = double(a) * double(b) + double(c); return float(s); }
+// Up to 2026-10-19 every fused operation went through f64 (the product exact there, the sum rounded to 53 bits, then
+// to f32), which this GPU runs at a small fraction of the f32 rate. They are the driver's f32 fma now, with two
+// exceptions that keep every result bit for bit what it was, on all 2^32 arguments (RX 9070 XT, RADV; NaN payloads
+// included):
+// - the polynomial stays in f64: the driver drops the f32 rounding between nested f64 operations (it folds
+//   f2f64(f2f32(x)) to x; the conversions are not NoContraction), so its four nested steps always ran in f64 with one
+//   rounding at the end. In f32 (true fma), 24,836 arguments round one f32 step differently, 2 of them in binary16.
+// - exp's error term fma(w, c, -ph) next to ph = w * c: the driver makes that the product it already has and a sum,
+//   exactly 0, even with every operand precise or opaque. It is Dekker's exact product error: w split into two halves
+//   of 12 bits (w * 4097 = w * (2^12 + 1)), log2(e) into 1.4423828125 + 0.0003122091293334961 (11 and 12 bits), so
+//   every partial product is exact and so is each sum.
+// A NaN argument keeps its payload, quieted and positive, as the f64 conversion left it.
 float tanh_ocml(float x)
 {
     const float y = abs(x);
     precise float z;
     if (y < 0.625) {
         precise float y2 = y * y;
-        precise float p = fmad(y2, fmad(y2, fmad(y2, fmad(y2, -0.005700020585209131, 0.02063407190144062), -0.053737930953502655), 0.13331416249275208), -0.3333328068256378);
+        precise double y2d = double(y2);
+        precise double pd = (((y2d * -0.005700020585209131 + 0.02063407190144062) * y2d - 0.053737930953502655) * y2d +
+                             0.13331416249275208) * y2d - 0.3333328068256378;
+        precise float p = float(pd);
         precise float yp = y * p;
-        z = fmad(y2, yp, y);
+        z = fma(y2, yp, y);
     } else {
         const float w = 2.0 * y;
         precise float ph = w * 1.4426950216293335;
-        precise float f0 = fmad(w, 1.4426950216293335, -ph);
-        precise float pl = fmad(w, 1.925962855864327e-08, f0);
+        precise float wt = w * 4097.0;
+        precise float wh = wt - (wt - w);
+        precise float wl = w - wh;
+        precise float f1 = wh * 1.4423828125 - ph;
+        precise float f2 = f1 + wh * 0.0003122091293334961;
+        precise float f3 = f2 + wl * 1.4423828125;
+        precise float f0 = f3 + wl * 0.0003122091293334961;
+        precise float pl = fma(w, 1.925962855864327e-08, f0);
         const float e = roundEven(ph);
         precise float d = ph - e;
         precise float a = d + pl;
         float t = ldexp(exp2(a), int(e));
         if (w > 88.72283172607422) t = uintBitsToFloat(0x7f800000u);
         precise float s1 = t + 1.0;
-        z = fmad(-2.0, rcp_a(s1), 1.0);
+        z = fma(-2.0, rcp_a(s1), 1.0);
+        if (isnan(z)) z = uintBitsToFloat(floatBitsToUint(y) | 0x00400000u);
     }
     return x < 0.0 || floatBitsToUint(x) == 0x80000000u ? -z : z;
 }

engine/runtime.h

@@ -273,7 +273,10 @@ struct RFrameSettings {
 };
 class RModel {
 public:
-    RModel(VkPhysicalDevice physical, VkDevice device, const std::string& directory, VkPipelineCache cache = VK_NULL_HANDLE);
+    // int64Atomics: the device has shaderBufferInt64Atomics enabled; preset F then measures the exposure in one pass
+    // (engine/rr/exposure_f.comp PASS 2) instead of two
+    RModel(VkPhysicalDevice physical, VkDevice device, const std::string& directory, VkPipelineCache cache = VK_NULL_HANDLE,
+           bool int64Atomics = false);
     ~RModel();
     RModel(const RModel&) = delete;
     RModel& operator=(const RModel&) = delete;

engine/runtime_m.inc

@@ -512,9 +512,12 @@ void MEngine::recordFrame(VkCommandBuffer cb, VkQueryPool stamps, uint32_t first
     auto& e = *impl; auto& m = e.model;
     if (e.external) throw std::runtime_error("this engine's network is external: record the front and the back");
     recordFront(cb, stamps, first);
+    // timing only: D4R_ENGINE_M_LAYER_STAMPS=1 also stamps after every network layer, from query first + 7 on
+    static const bool layerStamps = std::getenv("D4R_ENGINE_M_LAYER_STAMPS") != nullptr;
     for (unsigned l = 0; l < MLayers; ++l) {
         const uint32_t w = e.gw / 2 / MScale[l], h = e.gh / 2 / MScale[l];
         dispatch(cb, m.network[l], e.layerSets[l], (w + (MShift[l] ^ ((e.frame & 1u) * 2u)) + 7) / 8, (h + MShift[l] + 7) / 8); emptyStage(cb);
+        if (stamps && layerStamps) vkCmdWriteTimestamp(cb, VK_PIPELINE_STAGE_ALL_COMMANDS_BIT, stamps, first + 7 + l);
     }
     recordBack(cb, stamps, stamps ? first + 3 : 0);
 }
@@ -525,11 +528,24 @@ std::vector<MEngine::DebugResource> MEngine::debugResources() const {
     std::vector<DebugResource> r{{"enc0-params", e.enc0P.handle, VK_NULL_HANDLE, e.enc0P.size, 0, 0},
                                  {"tokens", e.tokens.handle, VK_NULL_HANDLE, e.tokens.size, 0, 0},
                                  {"head", e.dec1Out.handle, VK_NULL_HANDLE, e.dec1Out.size, 0, 0},
+                                 // the network's own activations, layer by layer
+                                 {"skip1", e.skip1.handle, VK_NULL_HANDLE, e.skip1.size, 0, 0},
+                                 {"merged1", e.merged1.handle, VK_NULL_HANDLE, e.merged1.size, 0, 0},
+                                 {"skip2", e.skip2.handle, VK_NULL_HANDLE, e.skip2.size, 0, 0},
+                                 {"merged2", e.merged2.handle, VK_NULL_HANDLE, e.merged2.size, 0, 0},
+                                 {"tubeA", e.tubeA.handle, VK_NULL_HANDLE, e.tubeA.size, 0, 0},
+                                 {"tubeB", e.tubeB.handle, VK_NULL_HANDLE, e.tubeB.size, 0, 0},
+                                 {"dec2", e.dec2Out.handle, VK_NULL_HANDLE, e.dec2Out.size, 0, 0},
                                  {"gates", VK_NULL_HANDLE, e.gates.handle, VkDeviceSize(w) * h * 4, w, h},
                                  {"mix", VK_NULL_HANDLE, e.mix.handle, VkDeviceSize(w) * h * 4, w, h},
                                  {"cov", VK_NULL_HANDLE, e.cov.handle, VkDeviceSize(w) * h * 8, w, h},
                                  {"mvA", VK_NULL_HANDLE, e.mvA[1 - e.current].handle, VkDeviceSize(w) * h * 4, w, h},
-                                 {"feat", VK_NULL_HANDLE, e.feat[1 - e.current].handle, VkDeviceSize(w) * h * 4, w, h}};
+                                 {"feat", VK_NULL_HANDLE, e.feat[1 - e.current].handle, VkDeviceSize(w) * h * 4, w, h},
+                                 // the reconstruction's outputs: the history, its half-resolution mean and the packed
+                                 // image the downsample reads
+                                 {"hist", VK_NULL_HANDLE, e.hist[1 - e.current].handle, VkDeviceSize(e.hw) * e.hh * 8, e.hw, e.hh},
+                                 {"half", VK_NULL_HANDLE, e.half[1 - e.current].handle, VkDeviceSize(e.hw / 2) * (e.hh / 2) * 16, e.hw / 2, e.hh / 2},
+                                 {"packed", e.packed.handle, VK_NULL_HANDLE, e.packed.size, 0, 0}};
     if (e.model.presetL) r.push_back({"skip", e.skipL.handle, VK_NULL_HANDLE, e.skipL.size, 0, 0});
     return r;
 }

engine/runtime_rr.inc

@@ -90,7 +90,7 @@ struct RModel::Impl {
     bool uploaded = false;
     Impl(VkPhysicalDevice physical, VkDevice dev) : device(physical, dev) {}
 };
-RModel::RModel(VkPhysicalDevice physical, VkDevice device, const std::string& directory, VkPipelineCache cache)
+RModel::RModel(VkPhysicalDevice physical, VkDevice device, const std::string& directory, VkPipelineCache cache, bool int64Atomics)
     : impl(std::make_unique<Impl>(physical, device)) {
     auto& m = *impl;
     // rr.txt: "D4RR0001", "preset N", one "stage NAME SHADER.spv SUBGROUP WEIGHTS TEXTURES SURFACES" line per stage in
@@ -103,7 +103,7 @@ RModel::RModel(VkPhysicalDevice physical, VkDevice device, const std::string& di
     if (line.rfind("D4RR0002", 0) == 0) {
         m.preset = 6;
         m.v2 = std::make_unique<R2Model>();
-        m.v2->load(m.device, directory, lines, cache);
+        m.v2->load(m.device, directory, lines, cache, int64Atomics);
         return;
     }
     if (line.rfind("D4RR0001", 0) != 0) throw std::runtime_error("not a Ray Reconstruction model (rr.txt D4RR0001)");

engine/runtime_rr2.inc

@@ -44,9 +44,10 @@ struct R2Model {
     Buffer aux;                   // encoder 0's second table (argument 104)
     std::vector<char> auxPending;
     std::array<Pipeline, 2> exposure;
+    Pipeline exposureOne;         // both exposure passes in one (exposure_f2.spv, 64-bit atomics); else null
     bool particles = false;       // the model has R2DownP
     bool gather = false;          // the model has R2Gather
-    void load(const Device& d, const std::string& directory, std::istringstream& lines, VkPipelineCache cache) {
+    void load(const Device& d, const std::string& directory, std::istringstream& lines, VkPipelineCache cache, bool int64Atomics) {
         const auto b = VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, t = VK_DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER, s = VK_DESCRIPTOR_TYPE_STORAGE_IMAGE;
         auto split = [](const std::string& list) {
             std::vector<std::string> items;
@@ -108,6 +109,8 @@ struct R2Model {
         blob(aux, auxPending, read(directory + "/enc0-aux.bin"));
         exposure[0].create(d.handle, directory + "/exposure_f0.spv", {b, t, b}, 64, cache);
         exposure[1].create(d.handle, directory + "/exposure_f1.spv", {b, b, b, t, s}, 64, cache);
+        if (int64Atomics && exists(directory + "/exposure_f2.spv"))
+            exposureOne.create(d.handle, directory + "/exposure_f2.spv", {b, t, b, b, t, s}, 64, cache);
     }
     void recordUpload(VkCommandBuffer cb) {   // the weights into the command buffer (recordContents)
         for (auto& st : stages) recordContents(cb, st.weights, st.pending);
@@ -172,7 +175,7 @@ struct R2Engine {
         std::vector<VkDescriptorSet> down;
         std::vector<std::array<VkDescriptorSet, 2>> downP;
         VkDescriptorSet pooled = VK_NULL_HANDLE, exposure0 = VK_NULL_HANDLE;
-        std::array<VkDescriptorSet, 2> exposure1{};
+        std::array<VkDescriptorSet, 2> exposure1{}, exposureOne{};
     };
     std::vector<Bindings> bindings;
     const Bindings* bound = nullptr;
@@ -337,7 +340,7 @@ struct R2Engine {
         exposure.blocksX = (rw + 15) / 16; exposure.blocks = exposure.blocksX * ((rh + 15) / 16);
         exposure.last[0] = static_cast<int32_t>(rw - 1); exposure.last[1] = static_cast<int32_t>(rh - 1);
         exposure.invCount = static_cast<float>(10000.0 / (4294967296.0 * rw * rh));
-        sums.create(d, exposure.blocks * 8, storage);
+        sums.create(d, std::max(exposure.blocks, 2u) * 8, storage);   // the one-pass exposure: the total and the count
         const size_t tiles = src.size();
         auto place = [&](size_t bytes) { const VkDeviceSize at = paramBytes; paramBytes += (bytes + 255) & ~VkDeviceSize(255); return at; };
         for (unsigned i = 0; i < R2Network; ++i) { blocks[i] = model.stages[i].args; paramAt[i] = place(blocks[i].size()); }
@@ -593,6 +596,17 @@ struct R2Engine {
                     b.downP.push_back({stageSet(b, R2DownP, downPAt[i], 0, output), stageSet(b, R2DownP, downPAt[i], 1, output)});
             }
             b.pooled = stageSet(b, R2Pooled, pooledAt, 0, output);
+            if (model.exposureOne.handle)
+                for (uint32_t par = 0; par < 2; ++par) {
+                    VkDescriptorSet set = b.exposureOne[par] = descriptorSet(d.handle, b.pool, model.exposureOne.descriptors);
+                    bindRange(d.handle, set, 0, paramAll, expAt);
+                    bindImage(d.handle, set, 1, images.color, nearest, images.layout);
+                    bindBuffer(d.handle, set, 2, sums);
+                    bindBuffer(d.handle, set, 3, values);
+                    bindImage(d.handle, set, 4, e[par ^ 1].view, nearest);
+                    bindImage(d.handle, set, 5, e[par].view);
+                }
+            else {
                 b.exposure0 = descriptorSet(d.handle, b.pool, model.exposure[0].descriptors);
                 bindRange(d.handle, b.exposure0, 0, paramAll, expAt);
                 bindImage(d.handle, b.exposure0, 1, images.color, nearest, images.layout);
@@ -605,6 +619,7 @@ struct R2Engine {
                     bindImage(d.handle, set, 3, e[par ^ 1].view, nearest);
                     bindImage(d.handle, set, 4, e[par].view);
                 }
+            }
             bindings.push_back(std::move(b));
             found = bindings.end() - 1;
         }
@@ -675,9 +690,15 @@ struct R2Engine {
             vkCmdUpdateBuffer(cb, paramAll.handle, at, std::min<VkDeviceSize>(65536, paramBytes - at), all.data() + at);
         memoryBarrier(cb, VK_PIPELINE_STAGE_TRANSFER_BIT, VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT, VK_ACCESS_TRANSFER_WRITE_BIT,
                       VK_ACCESS_SHADER_READ_BIT | VK_ACCESS_SHADER_WRITE_BIT);
+        // the exposure: one pass with 64-bit atomics (the last block finishes, as NGX's kernel), or the block sums and
+        // one workgroup that adds them
+        if (model.exposureOne.handle)
+            dispatch(cb, model.exposureOne, b.exposureOne[cur], exposure.blocksX, exposure.blocks / exposure.blocksX);
+        else {
             dispatch(cb, model.exposure[0], b.exposure0, exposure.blocksX, exposure.blocks / exposure.blocksX);
             step();
             dispatch(cb, model.exposure[1], b.exposure1[cur], 1, 1);
+        }
         step();
         stamp(1);
         for (unsigned stage = 0; stage < R2Network; ++stage) {

engine/test_m.py

@@ -90,7 +90,7 @@ def defines(name):
     # compile_m.py's arithmetic (ZLUDA's lowering); M_TEXL=0 gives the native accuracy kernels' for the network layers.
     # M_DEFS: extra defines for every layer_m.comp stage; M_DEFS_ENC1, M_DEFS_ENC3_TUBE, ...: for one layer
     texl = ['-DTEXL=1'] if kind >= 3 or os.environ.get('M_TEXL', '1') != '0' else []
-    # M_WPW=2: two 4x4-token windows per wave in the network layers (layer_m.comp WPW; not L's input / expansion)
+    # M_WPW=2: two 4x4-token windows per wave in the network layers (layer_m.comp WPW; L's input / expansion: M_DEFS_*)
     wpw = [f"-DWPW={os.environ['M_WPW']}"] if kind < 3 and os.environ.get('M_WPW') else []
     return [f'-DC={C}', f'-DNH={heads}', f'-DKIND={kind}', f'-DXC={NPM or CIN or 16}'] + texl + wpw + os.environ.get('M_DEFS', '').split() + os.environ.get('M_DEFS_' + name.upper(), '').split()
 

packaging/build_engine_models.sh

@@ -98,7 +98,7 @@ if (( DO_RR )); then
   done
   run ptx-f.log python3 "$KIT/kernels/tools/extract_dlss_ptx.py" "$DLL" "$WORK/ptx-f"
   run rr-f.log python3 "$KIT/engine/compile_rr2.py" --dll "$DLL" "$WORK/rr-f" --ptx-dir "$WORK/ptx-f"
-  place rr-f "$WORK/rr-f" $(model_files "$WORK/rr-f") enc0-aux.bin exposure_f0.spv exposure_f1.spv
+  place rr-f "$WORK/rr-f" $(model_files "$WORK/rr-f") enc0-aux.bin exposure_f0.spv exposure_f1.spv exposure_f2.spv
 fi
 
 if (( DO_SR )); then

tools/d4r_nvngx_shim.cpp

@@ -3385,6 +3385,9 @@ struct VulkanInterop
     // VK_EXT_shader_float8 (vkd3d-proton enables it when the GPU has it): M's and L's engine layers are FP8
     // cooperative-matrix shaders, which RADV's compiler crashes on without it (RDNA3)
     bool float8 = false;
+    // shaderBufferInt64Atomics is enabled on vkd3d-proton's device (it enables it when the GPU has it): Ray
+    // Reconstruction preset F then measures the exposure in one pass, as NGX does, instead of two
+    bool int64Atomics = false;
     // The engine runs on the game's own Vulkan device (native Vulkan NGX frontend): no vkd3d interop objects.
     bool native = false;
     PFN_vkGetPhysicalDeviceFormatProperties formatProperties = nullptr;
@@ -9322,6 +9325,17 @@ static void init_engine_vulkan()
     }
     logf("engine backend: VK_EXT_shader_float8 %s", g_vk.float8 ? "available (M and L engine layers)" : "unavailable: "
          "M and L run on the CUDA backend");
+    // the features vkd3d-proton created its device with
+    const VkPhysicalDeviceFeatures2* features = nullptr;
+    if (SUCCEEDED(g_vk.interop->GetDeviceFeatures(&features)) && features != nullptr)
+        for (auto* next = static_cast<const VkBaseInStructure*>(features->pNext); next != nullptr; next = next->pNext)
+        {
+            if (next->sType == VK_STRUCTURE_TYPE_PHYSICAL_DEVICE_VULKAN_1_2_FEATURES)
+                g_vk.int64Atomics = reinterpret_cast<const VkPhysicalDeviceVulkan12Features*>(next)->shaderBufferInt64Atomics == VK_TRUE;
+            if (next->sType == VK_STRUCTURE_TYPE_PHYSICAL_DEVICE_SHADER_ATOMIC_INT64_FEATURES)
+                g_vk.int64Atomics = reinterpret_cast<const VkPhysicalDeviceShaderAtomicInt64Features*>(next)->shaderBufferInt64Atomics == VK_TRUE;
+        }
+    logf("engine backend: 64-bit buffer atomics %s", g_vk.int64Atomics ? "enabled" : "not enabled: preset F's exposure in two passes");
 }
 
 static void init_engine_backend()
@@ -10316,7 +10330,7 @@ static void init_engine_backend_rr(int slot)
     const auto start = std::chrono::steady_clock::now();
     try
     {
-        auto model = std::make_unique<d4r::RModel>(g_vk.physical, g_vk.device, path);
+        auto model = std::make_unique<d4r::RModel>(g_vk.physical, g_vk.device, path, VK_NULL_HANDLE, g_vk.int64Atomics);
         if (model->preset() != g_rrPreset[slot])
             throw std::runtime_error("the directory holds a model of another preset");
         g_engineR.models[slot] = std::move(model);