verify window: PDL on sm_90+, mixer graph branches, batched Q4_0 query path, PLE post-ops in one pass (all opt-in)

Ported from #904 onto the 0.1.40 verify window. Every switch is OFF by default, so
default output and the default launch sequence are unchanged:

- STRATA_DF_PDL=1: programmatic dependent launch (pdl.hpp) for the window's
  quantizations and dense projections (multi-column MMVQ, multi-row BF16 GEMV);
  sm_90+ only, CUDA 12.3+; HIP and older cards keep ordinary launches.
- STRATA_DF_PLE=1: the window's hyper-connection write and PLE post-ops in one
  pass, which also writes the commit's history snapshots (native_gr_post_multi
  with strides, native_ple_postops_batch_snap).
- STRATA_DF_QB_Q4=1: a Q4_0 KV cache takes the batched QSA query path too.
- STRATA_DF_BRANCH=1: the mixer's side branches inside the window graph.
  Left opt-in because the author saw a stall on Linux with NVML queries.

The q/gate split now uses copy_rows_strided (a kernel) instead of a
cudaMemcpy2DAsync so a PDL chain stays kernel to kernel; the bytes are the same.

Left out: STRATA_DF_KVAPP (kv_append_q4_step_multi). 0.1.40 already appends a
window's K/V cells in one launch (kv_append_q4_steps, #783 PR-d).

Merged by hand with the 0.1.40 forms of the same code: the exact-T BF16
multi-token kernel, the PAIR/TS MMVQ variants, QFUSE (q8_attn), the fused
norm+rope and the batched residual write.

Co-authored-by: Francesco Albano <[email protected]>

(PR #904)

CMakeLists.txt

@@ -728,6 +728,15 @@ if(STRATA_ENABLE_CUDA OR STRATA_ENABLE_HIP)
     target_link_libraries(mmvq_il_parity PRIVATE strata_kernels ${_strata_gpu_runtime_target})
     add_test(NAME mmvq_il_parity COMMAND mmvq_il_parity)
   endif()
+  # programmatic dependent launch (pdl.hpp) in a captured verify-shaped chain: bitwise vs ordinary launches over
+  # replays with poisoned buffers, and only kernel-to-kernel programmatic edges (GPU, synthetic, no model; the edges
+  # are only made on sm_90+, elsewhere it checks that there are none).  NVIDIA only: PDL is a CUDA feature.
+  if(STRATA_ENABLE_CUDA AND NOT STRATA_HIP_GFX906)
+    add_executable(pdl_parity src/kernels/pdl_parity.cpp)
+    target_include_directories(pdl_parity PRIVATE ${_strata_gpu_include_directories})
+    target_link_libraries(pdl_parity PRIVATE strata_kernels ${_strata_gpu_runtime_target})
+    add_test(NAME pdl_parity COMMAND pdl_parity)
+  endif()
   # the IQ kernels that decode a weight part once for every column/entry, bitwise vs the per-column ones
   # (GPU, synthetic, no model; --bench adds old/new timings)
   add_executable(iq_multi_parity src/kernels/iq_multi_parity.cpp)
@@ -935,6 +944,12 @@ if(STRATA_ENABLE_CUDA OR STRATA_ENABLE_HIP)
   target_include_directories(kv_q4_parity PRIVATE ${_strata_gpu_include_directories})
   target_link_libraries(kv_q4_parity PRIVATE strata_kernels ${_strata_gpu_runtime_target})
   add_test(NAME kv_q4_parity COMMAND kv_q4_parity)
+  # the verify window's one-launch forms (hyper-connection write, PLE post-ops with the commit's
+  # history snapshots) vs the per-token launches they replace, bitwise (GPU, synthetic, no model)
+  add_executable(verify_batch_parity src/kernels/verify_batch_parity.cpp)
+  target_include_directories(verify_batch_parity PRIVATE ${_strata_gpu_include_directories})
+  target_link_libraries(verify_batch_parity PRIVATE strata_kernels ${_strata_gpu_runtime_target})
+  add_test(NAME verify_batch_parity COMMAND verify_batch_parity)
   add_executable(cvec_parity src/kernels/cvec_parity.cpp)
   target_include_directories(cvec_parity PRIVATE ${_strata_gpu_include_directories})
   target_link_libraries(cvec_parity PRIVATE strata_kernels ${_strata_gpu_runtime_target})

include/strata/core/verify.hpp

@@ -397,6 +397,12 @@ private:
     cudaStream_t cs_ = nullptr;
     cudaStream_t sh_cs_ = nullptr;
     cudaEvent_t ev_fork_ = nullptr, ev_join_ = nullptr;
+    // STRATA_DF_BRANCH: a layer's mixer work that reads only the layer's input, captured as parallel graph branches
+    // on these side streams (record_window); the same kernels on the same inputs, only their order is freer
+    bool df_branch_ = false;
+    cudaStream_t df_side_[2] = {};
+    cudaEvent_t df_fork_ = nullptr;
+    cudaEvent_t df_join_[2] = {};
     cudaGraphExec_t exec_[9] = {};
     cudaGraphExec_t exec_nr_[9] = {};   // #871: the doorbell variant of a stage that is all-resident otherwise
     cudaGraphExec_t commit_exec_ = nullptr;

include/strata/kernels/gr.hpp

@@ -37,6 +37,8 @@ void gr_set_fp32_activations(bool enabled);
 /// Layer reads use fused HC pre; the final mixer (null w_inject) retains its unfused graph order.
 /// Subsequent gr_write calls select pinned sigmoid/scatter arithmetic under the same setting.
 void gr_set_native_mmvf(bool enabled);
+/// Whether gr_write currently uses the native write (native_gr_post).
+bool gr_native_mmvf_enabled();
 
 /// Geometry of the hyper-connection.  The real model: n_embd = 2560, hc = 4, hc_lr = 320.
 struct GrShapes {

include/strata/kernels/native_gr_postops.hpp

@@ -27,6 +27,11 @@ void native_gr_pre_gated_multi(const float* xn, float* gate, float* mixed,
 // output may equal residual exactly; all other overlaps are invalid.
 void native_gr_post(const float* residual, const float* block_out, const float* inject,
                     float* output, int n_embd, int hc, void* stream);
+/// native_gr_post for n_tok tokens in one launch, token t's operands at t times the given strides (floats): the same
+/// expression per element.  output may equal residual exactly.
+void native_gr_post_multi(const float* residual, const float* block_out, const float* inject, float* output,
+                          int n_embd, int hc, int n_tok, long long r_stride, long long b_stride, long long i_stride,
+                          void* stream);
 
 void native_gr_post_multi(const float* residual, const float* block_out, const float* inject,
                           float* output, int n_embd, int hc, int n_tok, void* stream);

include/strata/kernels/native_ple_postops.hpp

@@ -35,4 +35,11 @@ void native_ple_postops(const float* projected_key, const float* hidden,
 void native_ple_postops_batch(float* key, float* hidden, const float* value, float* history, const PleWeights& w,
                               float* query_norm, float* gated, float* gate, int T, void* stream);
 
+/// The verify window: native_ple_postops_batch whose norms divide by the column count at run time, exactly as
+/// native_gr_rms_norm_weighted (the per-token path's) does, so every value is bitwise that of T native_ple_postops
+/// calls with the history advanced after each; it also writes `snap` (T x 9 x 10240): the history after each token,
+/// as the per-token path's copies leave it for the commit.
+void native_ple_postops_batch_snap(float* key, float* hidden, const float* value, float* history, const PleWeights& w,
+                                   float* query_norm, float* gated, float* gate, int T, float* snap, void* stream);
+
 } // namespace strata::kernels

include/strata/kernels/pdl.hpp

@@ -0,0 +1,90 @@
+// include/strata/kernels/pdl.hpp - programmatic dependent launch (PDL) for the verify window's kernel chain.
+//
+// On sm_90+ a kernel launched with cudaLaunchAttributeProgrammaticStreamSerialization may start while the kernel
+// before it is still running: everything before `pdl_wait()` overlaps that kernel (and its launch latency), everything
+// after it sees the kernel's writes, exactly as after an ordinary launch.  The kernels that opt in only LOAD constant
+// weights before `pdl_wait()` and write nothing, so every value they compute is unchanged.  `pdl_trigger()` lets the
+// next kernel launch early (it still waits for this grid to finish before it reads anything).  Their activations and
+// outputs are not `__restrict__` on CUDA (STRATA_PDL_RESTRICT): with it the compiler may load them before the wait
+// (the reason llama.cpp's PDL drops __restrict__ too); the weights keep it.
+//
+// Stream capture turns the attribute into a programmatic graph edge, which CUDA allows only between two kernel nodes.
+// `launch_pdl` therefore uses it only while `pdl_scope()` is set (the verify window's recorder, on a device where
+// `pdl_supported()`), only for a kernel whose code was built for sm_90+ (PTX of an older architecture, JIT-compiled,
+// has no griddepcontrol.wait: such a kernel would race), and only when every node the capturing stream's next node
+// depends on is a kernel - never after a memcpy, a memset, a host node or an event node.  After a join of branches
+// (several kernels) each edge is programmatic and the wait covers them all; STRATA_DF_PDL=2 allows one predecessor
+// only.  Everything else is an ordinary launch.  On HIP all of this is an ordinary launch, the device functions are
+// no-ops and the kernels keep their plain loops (kPdlPrefetch).
+#pragma once
+
+#include <cuda_runtime.h>
+
+#include <utility>
+
+namespace strata::kernels {
+
+/// Whether launches on this thread may carry the PDL attribute (set by the verify window's recorder).
+bool& pdl_scope();
+
+/// The current device runs PDL: a CUDA 12.3+ build, sm_90+ (STRATA_EMULATE_CC counts), code built for it, and
+/// STRATA_DF_PDL is set (1, or 2 for single-predecessor edges only).  Always false on HIP.
+bool pdl_supported();
+
+#if (defined(__CUDACC__) || defined(__HIPCC__))
+
+// not __CUDA_ARCH__-dependent: the host pass must see the same kernel signature (MSVC rejects the template stubs
+// otherwise), so on CUDA these parameters lose __restrict__ for every architecture
+#if defined(__HIPCC__)
+#define STRATA_PDL_RESTRICT __restrict__
+#else
+#define STRATA_PDL_RESTRICT
+#endif
+
+#if defined(__HIPCC__)
+inline constexpr bool kPdlPrefetch = false;   // HIP: no PDL, the kernels keep their plain loops
+#else
+inline constexpr bool kPdlPrefetch = true;    // CUDA: the weights a kernel can load before pdl_wait() are loaded there
+/// `kernel` may be launched with the PDL attribute into `stream` now (see the header comment).
+bool pdl_launch_ok(const void* kernel, cudaStream_t stream);
+#endif
+
+__device__ __forceinline__ void pdl_wait() {
+#if !defined(__HIPCC__) && defined(__CUDA_ARCH__) && (__CUDA_ARCH__ >= 900)
+    asm volatile("griddepcontrol.wait;" ::: "memory");
+#endif
+}
+
+__device__ __forceinline__ void pdl_trigger() {
+#if !defined(__HIPCC__) && defined(__CUDA_ARCH__) && (__CUDA_ARCH__ >= 900)
+    asm volatile("griddepcontrol.launch_dependents;" :::);
+#endif
+}
+
+/// `kernel<<<grid, block, smem, stream>>>(args...)`, with the PDL attribute when `pdl_launch_ok` says so.  A launch
+/// error stays the thread's last error for the caller's own check (cudaGetLastError) to report.
+template <typename... KArgs, typename... Args>
+inline void launch_pdl(void (*kernel)(KArgs...), dim3 grid, dim3 block, size_t smem, cudaStream_t stream,
+                       Args&&... args) {
+#if !defined(__HIPCC__)
+    if (pdl_scope() && pdl_launch_ok((const void*) kernel, stream)) {
+        cudaLaunchConfig_t cfg = {};
+        cfg.gridDim = grid;
+        cfg.blockDim = block;
+        cfg.dynamicSmemBytes = smem;
+        cfg.stream = stream;
+        cudaLaunchAttribute attr[1];
+        attr[0].id = cudaLaunchAttributeProgrammaticStreamSerialization;
+        attr[0].val.programmaticStreamSerializationAllowed = 1;
+        cfg.attrs = attr;
+        cfg.numAttrs = 1;
+        (void) cudaLaunchKernelEx(&cfg, kernel, std::forward<Args>(args)...);
+        return;
+    }
+#endif
+    kernel<<<grid, block, smem, stream>>>(std::forward<Args>(args)...);
+}
+
+#endif  // __CUDACC__ || __HIPCC__
+
+}  // namespace strata::kernels

include/strata/kernels/verify_kernels.hpp

@@ -41,8 +41,13 @@ void gdn_step_norm_multi(float* state, const float* h, int conv_channels, const
 /// Spin until *flag >= value (a mapped host flag).  The value is fixed at capture, so several rings can be
 /// outstanding at once (the split verify window keeps two).
 void wait_flag_ge(const uint32_t* flag, uint32_t value, void* stream);
-/// the GPU's %globaltimer (ns) into buf[i] (a one-thread kernel: the verify window's stage profiler).
+/// the GPU's %globaltimer (ns) into buf[i] (a one-thread kernel: the verify window's stage profiler).  Inside a PDL
+/// stretch (pdl.hpp) the stamp passes the early launch on, so a profiled window keeps the chain it measures.
 void gpu_stamp(unsigned long long* buf, int i, void* stream);
+/// dst[r][0, w) = src[r][0, w) for `rows` rows of source stride `src_w` floats (the query half of each q/gate head
+/// pair): the strided device-to-device cudaMemcpy2DAsync as a kernel, so the window's chain stays kernel to kernel.
+/// w and src_w multiples of 4 floats, both pointers 16-byte aligned.
+void copy_rows_strided(float* dst, const float* src, int64_t rows, int64_t w, int64_t src_w, void* stream);
 
 // ---- perf-review E-6: a layer whose routed experts are all in VRAM needs nothing from the host
 /// One group's plan, built on the device when every routed expert of its n*k entries is resident: the host pool's

src/core/verify.cpp

@@ -23,10 +23,13 @@
 #include "strata/kernels/kv_q4.hpp"
 #include "strata/kernels/kv_q8.hpp"
 #include "strata/kernels/native_mmvq.hpp"
+#include "strata/kernels/native_gr_postops.hpp"
+#include "strata/kernels/native_ple_postops.hpp"
 #include "strata/kernels/native_qsa.hpp"
 #include "strata/kernels/native_qsa_indexer.hpp"
 #include "strata/kernels/native_rope.hpp"
 #include "strata/kernels/ngram.hpp"
+#include "strata/kernels/pdl.hpp"
 #include "strata/kernels/ple.hpp"
 #include "strata/kernels/qsa.hpp"
 #include "strata/kernels/qsa_decode_attn.hpp"
@@ -360,6 +363,11 @@ Verifier::~Verifier() {
     if (prof_pin_) cudaFreeHost(prof_pin_);
     if (ev_fork_) cudaEventDestroy(ev_fork_);
     if (ev_join_) cudaEventDestroy(ev_join_);
+    for (cudaStream_t s : df_side_)
+        if (s) cudaStreamDestroy(s);
+    if (df_fork_) cudaEventDestroy(df_fork_);
+    for (cudaEvent_t e : df_join_)
+        if (e) cudaEventDestroy(e);
     if (arena_) cudaFree(arena_);
     void* hosts[] = {h_tok_, h_step_, h_pos_, h_commit_, h_ple_, h_out_, h_x_, h_ids_, h_w_, h_seq_, h_flag_, h_ymiss_,
                      h_flagA_, h_plan_, h_flagB_, h_plan_err_};
@@ -570,6 +578,19 @@ bool Verifier::init(const WeightTable& wt, const ModelGeometry& g, SessionState&
         err = "verify: event create failed";
         return false;
     }
+    {   // STRATA_DF_BRANCH=1 (opt-in): the mixer's side branches.  Off by default: on Linux (open modules, GSP) a
+        // stage's window with the branches stopped inside its graph when an NVML query (nvidia-smi) came between
+        // requests (2x RTX 3090 and 2x RTX 2080 Ti, #905); never seen on Windows.  0 / unset: one stream, as before
+        const char* v = std::getenv("STRATA_DF_BRANCH");
+        df_branch_ = v != nullptr && std::atoi(v) != 0;
+        if (df_branch_ && df_fork_ == nullptr) {
+            bool okb = cudaEventCreateWithFlags(&df_fork_, cudaEventDisableTiming) == cudaSuccess;
+            for (int i = 0; i < 2 && okb; ++i)
+                okb = cudaStreamCreateWithFlags(&df_side_[i], cudaStreamNonBlocking) == cudaSuccess &&
+                      cudaEventCreateWithFlags(&df_join_[i], cudaEventDisableTiming) == cudaSuccess;
+            if (!okb) { cudaGetLastError(); df_branch_ = false; }
+        }
+    }
     // Check whether 100% of experts across [lb_, le_) are resident in this stage's VRAM cache.
     // When true, every layer plans on device and writes directly into parts_ without any CPU doorbells,
     // wait_flag_ge spins, PCIe empty launches, or moe_hit_add copies.
@@ -668,15 +689,45 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
     };
     bool sg_gated_[2] = {false, false};   // per group: the combine applies the shared gate
     if (!batch_rec_) groups_[T] = G;
+    // Programmatic dependent launch (pdl.hpp; sm_90+, opt-in: STRATA_DF_PDL=1): the window's quantizations and dense
+    // projections may start while the kernel before them finishes, loading their weights before they wait for its
+    // output.  Only a launch whose predecessors in the graph are all kernels gets programmatic edges; every value is
+    // the same as without them.
+    struct PdlScope {
+        bool prev;
+        explicit PdlScope(bool on) : prev(pdl_scope()) { pdl_scope() = on; }
+        ~PdlScope() { pdl_scope() = prev; }
+    } pdl_window(G == 1 && pdl_supported());
+    // STRATA_DF_BRANCH: mixer work that reads only the layer's input runs on a side stream inside the graph.  fork(i)
+    // makes side i depend on everything captured on `cs` so far, join(i) makes `cs` depend on everything captured on
+    // side i: a branch never races with what came before its fork or after its join.  Not in a split or batch window,
+    // nor under the stage profiler (whose stamps are on `cs`).
+    const bool br = df_branch_ && G == 1 && !batch_rec_ && !prof_on_;
+    auto fork = [&](int i) {
+        cudaEventRecord(df_fork_, cs);
+        cudaStreamWaitEvent(df_side_[i], df_fork_, 0);
+    };
+    auto join = [&](int i) {
+        cudaEventRecord(df_join_[i], df_side_[i]);
+        cudaStreamWaitEvent(cs, df_join_[i], 0);
+    };
     const int64_t nQall = g.n_qsa_layers();
     // a batch window: row t is slot t, whose state lives in its own session
     auto slot_ss = [&](int t) -> SessionState& { return batch_rec_ ? *slots_[(size_t) brow_[t]] : ss; };
     const int hrow0 = batch_rec_ ? row_base_ : 0;   // a slot group's own hand-off rows
 
-    // ---- the window's inputs, from mapped staging
+    // ---- the window's inputs, from mapped staging.  With branches the steps and positions (first read by the first
+    // layer's mixer) cross PCIe on side 1 beside the embedding and the first hyper-connection read.
     copy_i32_from_mapped(tok_, m_tok_, T, cs);
-    copy_i32_from_mapped(step_, m_step_, (int64_t) T * kStepCount, cs);
-    copy_i32_from_mapped(pos_, m_pos_, (int64_t) MT * (NH + NKV + IQ), cs);
+    bool inputs_pending = false;
+    cudaStream_t in_s = cs;
+    if (br) {
+        fork(1);
+        in_s = df_side_[1];
+        inputs_pending = true;
+    }
+    copy_i32_from_mapped(step_, m_step_, (int64_t) T * kStepCount, in_s);
+    copy_i32_from_mapped(pos_, m_pos_, (int64_t) MT * (NH + NKV + IQ), in_s);
     // per-ROW positions of the K rows [t][NKV] and the indexer query rows [t][IQ] (for batched RoPE)
     const int32_t* pos_k = pos_ + MT * NH;
     const int32_t* pos_i = pos_ + MT * (NH + NKV);
@@ -781,9 +832,35 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                     return false;
                 }
             }
+            // STRATA_DF_PLE (opt-in: STRATA_DF_PLE=1): the window's post-projection steps in one pass, bitwise the per-token loop
+            // below (the norms divide at run time as its norms do; the history advances token by token), which also
+            // writes the commit's history snapshots.  Scratch: layer 0's expert rows, combined by now, and the
+            // hyper-connection read's `lo_`, which the read below rewrites.
+            static const bool ple_post_env = [] {
+                const char* e = std::getenv("STRATA_DF_PLE");
+                return e != nullptr && std::atoi(e) != 0;
+            }();
+            const bool ple_post = ple_batch_kv && ple_post_env && !batch_rec_;
             // #783 PR-g (stuchapin909): the window's residual writes in one launch (STRATA_NO_MULTI_GR=1: per token)
-            const bool multi_write = dec_batch && !g_no_multi_gr && n > 1;
+            const bool multi_write = !ple_post && dec_batch && !g_no_multi_gr && n > 1;
             if (multi_write) gr_write_multi(Rt(tb), bo_ + tb * N, inj2_ + tb * HC, gs, Rt(tb), n, cs);
+            if (ple_post) {
+                float* qn = parts_ + (size_t) tb * K * N;           // n x 10240 (K * N = 25600 per row)
+                float* gated = hit_out_ + (size_t) tb * K * N;      // n x 10240
+                float* pgate = lo_ + (size_t) tb * g.hc_lr;         // n x 4
+                try {
+                    if (gr_native_mmvf_enabled())   // gr_write's native write for the n rows in one launch
+                        native_gr_post_multi(Rt(tb), bo_ + tb * N, inj2_ + tb * HC, Rt(tb), (int) N, (int) HC, n, HC * N,
+                                             N, HC, cs);
+                    else
+                        for (int t = tb; t < te; ++t) gr_write(Rt(t), bo_ + t * N, inj2_ + t * HC, gs, Rt(t), cs);
+                    native_ple_postops_batch_snap(xn_ + (size_t) tb * HC * N, Rt(tb), z_ + (size_t) tb * N, ss.ple.hist,
+                                                  ss.ple.w, qn, gated, pgate, n, hist_snap_ + (size_t) tb * HS, cs);
+                } catch (const std::exception& e) {
+                    err = std::string("verify PLE (window): ") + e.what();
+                    return false;
+                }
+            } else
             for (int t = tb; t < te; ++t) {
                 if (batched) { pw.pre_key = ple_key_ + (size_t) t * NG_HC_DIM; pw.pre_value = ple_val_ + t * N; }
                 if (!multi_write) gr_write(Rt(t), bo_ + t * N, inj2_ + t * HC, gs, Rt(t), cs);
@@ -829,6 +906,10 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                                 (int) (l * kProfPer + (half == 0 ? 27 : 30)));
         };
         const bool q8_attn = gr_read_group(0, pending, inj2_, inj_);   // true: xq_ holds mixed's q8_1 (STRATA_QFUSE)
+        if (inputs_pending) {   // the window's steps and positions, before the first mixer reads them
+            join(1);
+            inputs_pending = false;
+        }
         stamp(l, 1, grp);
         float* xm = mixed_ + tb * N;
         bool il_ready = false;   // xil_ holds xq_'s interleaved copy (reset whenever xq_ is rewritten)
@@ -863,6 +944,15 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                 float* gate = gate_L_ + (size_t) gi * MT * HV;
                 float* beta = beta_L_ + (size_t) gi * MT * HV;
                 if (!q8_attn) native_quantize_q8_1(xm, xq_, (int) N, n, cs);
+                if (br) {   // STRATA_DF_BRANCH: a/b and z beside q/k/v and the conv; all of them read only this layer's input
+                    fork(0);
+                    gdn_ab_multi(xm, (const uint16_t*) wa->data, (const uint16_t*) wb->data, (const float*) wdt->data,
+                                 (const float*) wsa->data, gate + (size_t) tb * HV, beta + (size_t) tb * HV, (int) N,
+                                 (int) HV, n, df_side_[0]);
+                    fork(1);
+                    native_mmvq(wg->native_type, wg->native_data, xq_, z_ + (size_t) tb * ZV, (int) N, (int) ZV, n,
+                                df_side_[1]);
+                }
                 mm(wqkv, qkv + (size_t) tb * C, (int) N, (int) C);
                 stamp(l, 2, grp);
                 if (batch_rec_) {   // contiguous rows may be proposals for the same slot
@@ -880,12 +970,17 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                 gdn_conv_l2_multi(conv, qkv, (const float*) wc->data, hb, (int) C, (int) (2 * HK), EPS, n, cs, tb,
                                   self_commit);
                 stamp(l, 3, grp);
+                if (br) {   // gate, beta and z, before the recurrence reads them
+                    join(0);
+                    join(1);
+                } else {
                 gdn_ab_multi(xm, (const uint16_t*) wa->data, (const uint16_t*) wb->data, (const float*) wdt->data,
                              (const float*) wsa->data, gate + (size_t) tb * HV, beta + (size_t) tb * HV, (int) N, (int) HV,
                              n, cs);
                 stamp(l, 4, grp);
                 mm(wg, z_ + (size_t) tb * ZV, (int) N, (int) ZV);
                 stamp(l, 5, grp);
+                }
                 // the recurrence from the untouched state over tokens [0, te); outputs only for this group's
                 if (batch_rec_) {   // each slot's recurrence over its own proposed-token group
                     for (int t = tb; t < te;) {
@@ -920,6 +1015,13 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                 if (!native_of(wq, v.name("attn_q.weight"), err) || !native_of(wk, v.name("attn_k.weight"), err) ||
                     !native_of(wv, v.name("attn_v.weight"), err) || !native_of(wo, v.name("attn_output.weight"), err))
                     return false;
+                auto norm_rope_on = [&](float* data, const WeightRef* norm, int rows, int cols, const int32_t* pos,
+                                        cudaStream_t sx) {
+                    if (native_qsa_enabled()) native_qsa_rms_norm_weighted(data, (const float*) norm->data, data, cols, rows, EPS, sx);
+                    else rms_norm_weighted(data, (const float*) norm->data, rows, cols, EPS, sx);
+                    if (native_rope_enabled()) native_rope_apply(data, data, rows, cols, (int) s.n_rot, rope_scaling(), pos, sx);
+                    else rope_neox_apply(data, data, rows, cols, (int) s.n_rot, st.cos_tab, st.sin_tab, pos, sx);
+                };
                 // #783 PR-f (stuchapin909): the per-head RMSNorm and the rope in one launch, bit-identical to the pair
                 // (rope_parity check 6); STRATA_NO_NORM_ROPE=1 keeps the two; off on HIP until its parity check passes
                 const bool fuse_nr = native_qsa_enabled() && native_rope_enabled() && native_norm_rope_usable((int) HD, (int) s.n_rot);
@@ -929,16 +1031,35 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                                                  rope_scaling(), pos, cs);
                         return;
                     }
-                    if (native_qsa_enabled()) native_qsa_rms_norm_weighted(data, (const float*) norm->data, data, cols, rows, EPS, cs);
-                    else rms_norm_weighted(data, (const float*) norm->data, rows, cols, EPS, cs);
-                    if (native_rope_enabled()) native_rope_apply(data, data, rows, cols, (int) s.n_rot, rope_scaling(), pos, cs);
-                    else rope_neox_apply(data, data, rows, cols, (int) s.n_rot, st.cos_tab, st.sin_tab, pos, cs);
+                    norm_rope_on(data, norm, rows, cols, pos, cs);
                 };
                 float* idx_raw = idx_raw_L_ + (size_t) qi * MT * ID;
                 // the per-token GEMVs / norms / RoPEs / copies of this layer as one launch over the
                 // window's rows each - row-wise identical arithmetic (STRATA_DEC_BATCH=0: token by token)
-                const bool qb = dec_batch && n > 1 && native_qsa_enabled() && native_rope_enabled() && !st.kv_q4;
+                // STRATA_DF_QB_Q4 (opt-in: STRATA_DF_QB_Q4=1): a Q4_0 KV cache too - its one per-token extra, the query's rotation, is
+                // row-wise and the batched path rotates (270650e); 0: Q4_0 token by token, as before
+                static const bool qb_q4 = [] {
+                    const char* e = std::getenv("STRATA_DF_QB_Q4");
+                    return e != nullptr && std::atoi(e) != 0;
+                }();
+                const bool qb = dec_batch && n > 1 && native_qsa_enabled() && native_rope_enabled() && (!st.kv_q4 || qb_q4);
                 if (!q8_attn) native_quantize_q8_1(xm, xq_, (int) N, n, cs);
+                // STRATA_DF_BRANCH: the query side (q and its norm/RoPE/rotation on side 0, the indexer query on side 1)
+                // beside the K/V side; the indexer query joins before the block scores, the query before attention
+                const bool qbr = br && qb;
+                if (qbr) {
+                    fork(0);
+                    native_mmvq(wq->native_type, wq->native_data, xq_, qfull_ + tb * NH * 2 * HD, (int) N,
+                                (int) (NH * 2 * HD), n, df_side_[0]);
+                    copy_rows_strided(qcur_ + tb * NH * HD, qfull_ + tb * NH * 2 * HD, (int64_t) n * NH, HD, 2 * HD,
+                                      df_side_[0]);
+                    norm_rope_on(qcur_ + tb * NH * HD, wqn, (int) (n * NH), (int) HD, pos_ + tb * NH, df_side_[0]);
+                    if (st.kv_rot) fwht256_inplace_cuda(qcur_ + tb * NH * HD, (int64_t) n * NH, df_side_[0]);   // <Hq, Hk> = <q, k>
+                    fork(1);
+                    bf16_gemv_fp32_mmvf_multi(mixed_ + tb * N, N, (const uint16_t*) wiq->data, qidx_ + tb * IQ * ID,
+                                              IQ * ID, N, IQ * ID, n, df_side_[1]);
+                    norm_rope_on(qidx_ + tb * IQ * ID, wiqn, (int) (n * IQ), (int) ID, pos_i + tb * IQ, df_side_[1]);
+                }
                 if (qb) bf16_gemv_fp32_mmvf_multi(mixed_ + tb * N, N, (const uint16_t*) wik->data, idx_raw + tb * ID, ID, N, ID, n, cs);
                 else for (int t = tb; t < te; ++t)
                     bf16_gemv_fp32_mmvf(mixed_ + t * N, (const uint16_t*) wik->data, idx_raw + t * ID, (int) N, (int) ID, cs);
@@ -1023,6 +1144,9 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                     }
                 }
                 stamp(l, 9, grp);
+                if (qbr) {
+                    join(1);   // the indexer query, for the block scores
+                } else {
                 mm(wq, qfull_ + tb * NH * 2 * HD, (int) N, (int) (NH * 2 * HD));
                 if (qb) {
                     if (fuse_nr) {   // the q/gate split reads q straight out of the q|gate rows (stride 2 * HD)
@@ -1030,11 +1154,8 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                                                  qcur_ + tb * NH * HD, (int) (n * NH), (int) HD, (int) s.n_rot, EPS,
                                                  rope_scaling(), pos_ + tb * NH, cs);
                     } else {
-                    if (cudaMemcpy2DAsync(qcur_ + tb * NH * HD, (size_t) HD * 4, qfull_ + tb * NH * 2 * HD, (size_t) HD * 2 * 4,
-                                          (size_t) HD * 4, (size_t) (n * NH), cudaMemcpyDeviceToDevice, cs) != cudaSuccess) {
-                        err = "verify: the q/gate split failed";
-                        return false;
-                    }
+                    // the q/gate split as a copy kernel: the window's chain stays kernel to kernel (PDL)
+                    copy_rows_strided(qcur_ + tb * NH * HD, qfull_ + tb * NH * 2 * HD, (int64_t) n * NH, HD, 2 * HD, cs);
                     norm_rope(qcur_ + tb * NH * HD, wqn, (int) (n * NH), (int) HD, pos_ + tb * NH);
                     }
                     if (st.kv_rot) fwht256_inplace_cuda(qcur_ + tb * NH * HD, (int64_t) n * NH, cs);   // <Hq, Hk> = <q, k>
@@ -1063,6 +1184,7 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                     norm_rope(qx, wiqn, (int) IQ, (int) ID, pos_ + t * NH);
                 }
                 }
+                }
                 stamp(l, 10, grp);
                 if (batch_rec_) {   // each row selects and attends over its own slot's K/V
                     for (int t = tb; t < te; ++t) {
@@ -1086,6 +1208,7 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                 // KV streaming: the n selections' blocks resident (device-side, inside the graph)
                 qsa_kv_resolve(st, *g_, sel_ + (size_t) tb * cap_, step_ + tb * kStepCount, n, cap_, cs);
                 stamp(l, 12, grp);
+                if (qbr) join(0);   // the query, for attention
                 const QsaAttnPools pools = qsa_attn_pools(st);
                 qsa_decode_attn_batch(qcur_ + tb * NH * HD, pools, sel_ + (size_t) tb * cap_, step_ + tb * kStepCount, cap_,
                                       s, attn_scratch_ + (size_t) tb * attn_scratch_floats_, attn_ + tb * NH * HD, n, cs);
@@ -1563,6 +1686,33 @@ bool Verifier::capture(int T, std::string& err) {
         std::fprintf(stderr, "strata verify: the %d-token window graph has %zu nodes:", T, nn);
         for (size_t i = 0; i < v.size() && i < 40; ++i) std::fprintf(stderr, " %d x %.60s;", v[i].first, v[i].second.c_str());
         std::fprintf(stderr, "\n");
+#if CUDART_VERSION >= 12030   // the programmatic (PDL) edges: CUDA allows them only from a kernel node to a kernel node
+        size_t ne = 0;
+#if CUDART_VERSION >= 13000
+        cudaGraphGetEdges(graph, nullptr, nullptr, nullptr, &ne);
+#else
+        cudaGraphGetEdges_v2(graph, nullptr, nullptr, nullptr, &ne);
+#endif
+        std::vector<cudaGraphNode_t> from(ne), to(ne);
+        std::vector<cudaGraphEdgeData> ed(ne);
+#if CUDART_VERSION >= 13000
+        const cudaError_t ge = cudaGraphGetEdges(graph, from.data(), to.data(), ed.data(), &ne);
+#else
+        const cudaError_t ge = cudaGraphGetEdges_v2(graph, from.data(), to.data(), ed.data(), &ne);
+#endif
+        size_t prog = 0, bad = 0;
+        for (size_t i = 0; ge == cudaSuccess && i < ne; ++i) {
+            if (ed[i].type != cudaGraphDependencyTypeProgrammatic) continue;
+            ++prog;
+            cudaGraphNodeType ta, tb2;
+            if (cudaGraphNodeGetType(from[i], &ta) != cudaSuccess || cudaGraphNodeGetType(to[i], &tb2) != cudaSuccess ||
+                ta != cudaGraphNodeTypeKernel || tb2 != cudaGraphNodeTypeKernel)
+                ++bad;
+        }
+        if (ge != cudaSuccess) cudaGetLastError();
+        std::fprintf(stderr, "strata verify: %zu edges, %zu programmatic (PDL), %zu of them not kernel to kernel%s\n", ne,
+                     prog, bad, ge != cudaSuccess ? " (the edge query failed)" : "");
+#endif
     }
 #endif
     const bool made = instantiate_evicting(exec_t, graph, {}, "verify: instantiate: ", err);

src/kernels/cuda/gr.cu

@@ -300,6 +300,7 @@ __global__ void gr_write_kernel(const float* __restrict__ R, const float* __rest
 
 void gr_set_fp32_activations(bool enabled) { fp32_activations = enabled; }
 void gr_set_native_mmvf(bool enabled) { native_mmvf = enabled; }
+bool gr_native_mmvf_enabled() { return native_mmvf; }
 
 size_t gr_workspace_init(const GrShapes& s, void* base, GrWorkspace& out) {
     const size_t hc_dim = (size_t) s.hc * (size_t) s.n_embd;

src/kernels/cuda/native_bf16.cu

@@ -1,6 +1,7 @@
 #include "strata/kernels/bf16_gemv.hpp"
 #include "s26_tsum.cuh"
 #include "strata/kernels/bf16_bits.hpp"
+#include "strata/kernels/pdl.hpp"
 
 #include <cuda_runtime.h>
 #include <cstdlib>
@@ -74,8 +75,8 @@ __global__ void bf16_f32_mmvf_kernel(const float* __restrict__ x, const uint16_t
 // token keeps its own accumulator with exactly the single-row kernel's order (pairs, two ordered FMAs, the same warp
 // and block reductions), so each output is bit-identical to a bf16_f32_mmvf_kernel launch of its own.
 template <int BLOCK_SIZE, int NT, bool EXACT_T = true>
-__global__ void bf16_f32_mmvf_multi_kernel(const float* __restrict__ x, int64_t ldx, const uint16_t* __restrict__ w,
-                                          float* __restrict__ y, int64_t ldy, int n_in, int n_tok) {
+__global__ void bf16_f32_mmvf_multi_kernel(const float* STRATA_PDL_RESTRICT x, int64_t ldx, const uint16_t* __restrict__ w,
+                                          float* STRATA_PDL_RESTRICT y, int64_t ldy, int n_in, int n_tok) {
     const int t = threadIdx.x;
     const uint16_t* row = w + (size_t) blockIdx.x * n_in;
     const uint32_t* weights2 = reinterpret_cast<const uint32_t*>(row);
@@ -88,8 +89,18 @@ __global__ void bf16_f32_mmvf_multi_kernel(const float* __restrict__ x, int64_t
     float acc[NT];
 #pragma unroll
     for (int k = 0; k < NT; ++k) acc[k] = 0.0f;
-    for (int pair = t; pair < n_in / 2; pair += BLOCK_SIZE) {
-        const uint32_t weight = __ldg(weights2 + pair);
+    // PDL (pdl.hpp): this thread's first PRE weight pairs are loaded before waiting for the activations
+    constexpr int PRE = kPdlPrefetch ? 8 : 0;
+    uint32_t wpre[PRE > 0 ? PRE : 1];
+#pragma unroll
+    for (int q = 0; q < PRE; ++q) {
+        const int pair = t + q * BLOCK_SIZE;
+        wpre[q] = pair < n_in / 2 ? __ldg(weights2 + pair) : 0u;
+    }
+    pdl_wait();
+    int q = 0;
+    for (int pair = t; pair < n_in / 2; pair += BLOCK_SIZE, ++q) {
+        const uint32_t weight = q < PRE ? wpre[q] : __ldg(weights2 + pair);
         const float w0 = f32_from_bf16((uint16_t) weight), w1 = f32_from_bf16((uint16_t) (weight >> 16));
 #pragma unroll
         for (int k = 0; k < NT; ++k) {
@@ -306,13 +317,13 @@ void bf16_gemv_fp32_mmvf_multi(const float* x, int64_t ldx, const uint16_t* w, f
     }
 #define STRATA_MMVF_M(N) case N: \
     switch (n_tok) { \
-        case 1: bf16_f32_mmvf_multi_kernel<N, 1, true><<<(unsigned) n_out, N, 0, st>>>(x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
-        case 2: bf16_f32_mmvf_multi_kernel<N, 2, true><<<(unsigned) n_out, N, 0, st>>>(x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
-        case 3: bf16_f32_mmvf_multi_kernel<N, 3, true><<<(unsigned) n_out, N, 0, st>>>(x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
-        case 4: bf16_f32_mmvf_multi_kernel<N, 4, true><<<(unsigned) n_out, N, 0, st>>>(x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
-        case 5: bf16_f32_mmvf_multi_kernel<N, 5, true><<<(unsigned) n_out, N, 0, st>>>(x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
-        case 6: bf16_f32_mmvf_multi_kernel<N, 6, true><<<(unsigned) n_out, N, 0, st>>>(x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
-        default: bf16_f32_mmvf_multi_kernel<N, 8, false><<<(unsigned) n_out, N, 0, st>>>(x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
+        case 1: launch_pdl(bf16_f32_mmvf_multi_kernel<N, 1, true>, dim3((unsigned) n_out), dim3(N), 0, st, x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
+        case 2: launch_pdl(bf16_f32_mmvf_multi_kernel<N, 2, true>, dim3((unsigned) n_out), dim3(N), 0, st, x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
+        case 3: launch_pdl(bf16_f32_mmvf_multi_kernel<N, 3, true>, dim3((unsigned) n_out), dim3(N), 0, st, x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
+        case 4: launch_pdl(bf16_f32_mmvf_multi_kernel<N, 4, true>, dim3((unsigned) n_out), dim3(N), 0, st, x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
+        case 5: launch_pdl(bf16_f32_mmvf_multi_kernel<N, 5, true>, dim3((unsigned) n_out), dim3(N), 0, st, x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
+        case 6: launch_pdl(bf16_f32_mmvf_multi_kernel<N, 6, true>, dim3((unsigned) n_out), dim3(N), 0, st, x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
+        default: launch_pdl(bf16_f32_mmvf_multi_kernel<N, 8, false>, dim3((unsigned) n_out), dim3(N), 0, st, x, ldx, w, y, ldy, (int) n_in, n_tok); break; \
     } break
     switch (mmvf_block_size(n_in)) {
         STRATA_MMVF_M(32); STRATA_MMVF_M(64); STRATA_MMVF_M(96); STRATA_MMVF_M(128);

src/kernels/cuda/native_gr_postops.cu

@@ -81,6 +81,19 @@ __global__ void post(const float* residual, const float* __restrict__ block_out,
     // Exact residual/output alias is supported; no other thread reads residual[i].
     output[i] = __fmaf_rn(block_out[d], weight, residual[i]);
 }
+// `post` for n tokens (blockIdx.y = token): the same expression per element
+__global__ void post_multi(const float* residual, const float* __restrict__ block_out,
+                           const float* __restrict__ inject, float* output, int n_embd, int hc, float scale,
+                           long long r_stride, long long b_stride, long long i_stride) {
+    const std::size_t i = std::size_t(blockIdx.x) * blockDim.x + threadIdx.x;
+    if (i >= std::size_t(n_embd) * hc) return;
+    const long long t = blockIdx.y;
+    residual += t * r_stride; output += t * r_stride; block_out += t * b_stride; inject += t * i_stride;
+    const int c = int(i / n_embd), d = int(i % n_embd);
+    const float weight = scale_zero_bias(sigmoid(scale_zero_bias(inject[c], scale)), 2.0f);
+    // Exact residual/output alias is supported; no other thread reads residual[i].
+    output[i] = __fmaf_rn(block_out[d], weight, residual[i]);
+}
 void check_pointer(const void* p) {
     if (!p || reinterpret_cast<std::uintptr_t>(p) % alignof(float))
         throw std::invalid_argument("native GR postops require non-null four-byte aligned pointers");
@@ -135,4 +148,14 @@ void native_gr_post_multi(const float* residual, const float* block_out, const f
         residual, block_out, inject, output, n_embd, hc, 1.0f / float(hc));
     check_launch();
 }
+void native_gr_post_multi(const float* residual, const float* block_out, const float* inject, float* output,
+                          int n_embd, int hc, int n_tok, long long r_stride, long long b_stride, long long i_stride,
+                          void* stream) {
+    check_shape(n_embd, hc);
+    check_pointer(residual); check_pointer(block_out); check_pointer(inject); check_pointer(output);
+    if (n_tok < 1) return;
+    post_multi<<<dim3(blocks(std::size_t(n_embd) * hc), unsigned(n_tok)), THREADS, 0, static_cast<cudaStream_t>(stream)>>>(
+        residual, block_out, inject, output, n_embd, hc, 1.0f / float(hc), r_stride, b_stride, i_stride);
+    check_launch();
+}
 } // namespace strata::kernels

src/kernels/cuda/native_mmvq.cu

@@ -31,6 +31,7 @@
 #if !defined(__HIPCC__)
 #include "q8_1_il.cuh"
 #endif
+#include "strata/kernels/pdl.hpp"
 
 #include <cuda_fp16.h>
 #include <cuda_runtime.h>
@@ -149,8 +150,10 @@ __device__ __forceinline__ float warp_max(float x) {
 }
 
 __launch_bounds__(QUANT_THREADS, 1)
-__global__ void native_quantize_q8_1_kernel(const float* __restrict__ x,
-                                           Q81Block* __restrict__ y, int n_in) {
+__global__ void native_quantize_q8_1_kernel(const float* STRATA_PDL_RESTRICT x,
+                                           Q81Block* STRATA_PDL_RESTRICT y, int n_in) {
+    pdl_trigger();   // PDL (pdl.hpp): the projection after this one may start loading its weights
+    pdl_wait();
     const int i = int(blockIdx.x) * QUANT_THREADS + int(threadIdx.x);
     if (i >= n_in) return; // n_in is a multiple of 32: only whole warps return.
     const float xi = x[i];
@@ -1160,8 +1163,8 @@ static bool s26_tsum_on() {
 template<typename F, int NCOLS, int NW, int ROWS, bool TS = false, bool PAIR = false>
 __launch_bounds__(NW * WARP, (ROWS <= 2 ? 4 : 1))
 __global__ void native_mmvq_multi_kernel(const typename F::Block* __restrict__ w,
-                                         const Q81Block* __restrict__ x,
-                                         float* __restrict__ y, int n_in, int n_out,
+                                         const Q81Block* STRATA_PDL_RESTRICT x,
+                                         float* STRATA_PDL_RESTRICT y, int n_in, int n_out,
                                          const typename F::Block* __restrict__ w2 = nullptr, float* __restrict__ y2 = nullptr) {
     constexpr int BPI = F::BPI * NW / WARPS;           // blocks per iteration scale with the warp count
     const int tid = WARP * int(threadIdx.y) + int(threadIdx.x);
@@ -1173,7 +1176,33 @@ __global__ void native_mmvq_multi_kernel(const typename F::Block* __restrict__ w
     const int blocks_per_row = n_in / F::DIV;
     const int x_stride = n_in / Q8K;                   // Q8_1 blocks per activation column
     float tmp[NCOLS][ROWS] = {};
-    for (int kbx = tid / F::T; kbx < blocks_per_row; kbx += BPI) {
+    int kbx = tid / F::T;
+    if constexpr (kPdlPrefetch) {
+        // PDL (pdl.hpp): the first block's weights are loaded before waiting for the activations, then that block
+        // is applied exactly as the loop below would apply it (the same rows, columns and order)
+        typename F::W w0[ROWS];
+        const int kqs = F::kqs(tid);
+        if (kbx < blocks_per_row) {
+#pragma unroll
+            for (int i = 0; i < ROWS; ++i)
+                if (row0 + i < n_out) w0[i] = F::load(w + std::size_t(row0 + i) * blocks_per_row + kbx, kqs);
+        }
+        pdl_trigger();
+        pdl_wait();
+        if (kbx < blocks_per_row) {
+            const int kby = kbx * F::KBY;
+#pragma unroll
+            for (int i = 0; i < ROWS; ++i) {
+                if (row0 + i < n_out) {
+#pragma unroll
+                    for (int j = 0; j < NCOLS; ++j)
+                        tmp[j][i] += F::apply(w0[i], x + std::size_t(j) * x_stride + kby, kqs);
+                }
+            }
+        }
+        kbx += BPI;
+    }
+    for (; kbx < blocks_per_row; kbx += BPI) {
         const int kby = kbx * F::KBY;
         const int kqs = F::kqs(tid);
 #pragma unroll
@@ -1234,10 +1263,14 @@ void launch_multi_n(const void* weights, const void* x_q8_1, float* y, int n_in,
     if (!g_multi_exact) {
         constexpr int NW = NCOLS <= 4 ? 4 : 2;
         const unsigned blocks = unsigned((std::size_t(n_out) + 1) / 2);
-        native_mmvq_multi_kernel<F, NCOLS, NW, 2><<<blocks, dim3(WARP, NW), 0, s>>>(w, x, y, n_in, n_out);
+        launch_pdl(native_mmvq_multi_kernel<F, NCOLS, NW, 2, false, false>, dim3(blocks), dim3(WARP, NW), 0, s, w, x, y, n_in, n_out,
+                   (const typename F::Block*) nullptr, (float*) nullptr);
         return;
     }
     const dim3 threads(WARP, WARPS);
+    // the PAIR arguments, unused here (launch_pdl passes every parameter)
+    const typename F::Block* const no_w2 = nullptr;
+    float* const no_y2 = nullptr;
     // #783 PR-h (stuchapin909): two rows per block for every K, not only small K - each row keeps its own partial sums and
     // reduction, so a row's result does not depend on its neighbour; STRATA_NO_MMVQ_ROWS2=1 keeps one row for large K
     static const bool rows1 = [] {
@@ -1245,14 +1278,14 @@ void launch_multi_n(const void* weights, const void* x_q8_1, float* y, int n_in,
         return v != nullptr && v[0] != '\0' && v[0] != '0';
     }();
     if (rows1 && n_in / F::DIV >= F::BPI) {
-        if (s26_tsum_on()) native_mmvq_multi_kernel<F, NCOLS, WARPS, 1, true><<<unsigned(n_out), threads, 0, s>>>(w, x, y, n_in, n_out);
-        else native_mmvq_multi_kernel<F, NCOLS, WARPS, 1><<<unsigned(n_out), threads, 0, s>>>(w, x, y, n_in, n_out);
+        if (s26_tsum_on()) launch_pdl(native_mmvq_multi_kernel<F, NCOLS, WARPS, 1, true, false>, dim3(unsigned(n_out)), threads, 0, s, w, x, y, n_in, n_out, no_w2, no_y2);
+        else launch_pdl(native_mmvq_multi_kernel<F, NCOLS, WARPS, 1, false, false>, dim3(unsigned(n_out)), threads, 0, s, w, x, y, n_in, n_out, no_w2, no_y2);
         return;
     }
     constexpr int ROWS = 2;
     const unsigned blocks = unsigned((std::size_t(n_out) + ROWS - 1) / ROWS);
-    if (s26_tsum_on()) native_mmvq_multi_kernel<F, NCOLS, WARPS, ROWS, true><<<blocks, threads, 0, s>>>(w, x, y, n_in, n_out);
-    else native_mmvq_multi_kernel<F, NCOLS, WARPS, ROWS><<<blocks, threads, 0, s>>>(w, x, y, n_in, n_out);
+    if (s26_tsum_on()) launch_pdl(native_mmvq_multi_kernel<F, NCOLS, WARPS, ROWS, true, false>, dim3(blocks), threads, 0, s, w, x, y, n_in, n_out, no_w2, no_y2);
+    else launch_pdl(native_mmvq_multi_kernel<F, NCOLS, WARPS, ROWS, false, false>, dim3(blocks), threads, 0, s, w, x, y, n_in, n_out, no_w2, no_y2);
 }
 
 template<typename F>
@@ -2042,8 +2075,8 @@ void native_quantize_q8_1(const float* x, void* x_q8_1, int n_in, int ncols, voi
     // ncols * n_in elements: every 32-element block stays inside one column.
     const int n_total = n_in * ncols;
     const unsigned blocks = unsigned((std::size_t(n_total) + QUANT_THREADS - 1) / QUANT_THREADS);
-    native_quantize_q8_1_kernel<<<blocks, QUANT_THREADS, 0,
-                                 static_cast<cudaStream_t>(stream)>>>(x, static_cast<Q81Block*>(x_q8_1), n_total);
+    launch_pdl(native_quantize_q8_1_kernel, dim3(blocks), dim3(QUANT_THREADS), 0, static_cast<cudaStream_t>(stream), x,
+               static_cast<Q81Block*>(x_q8_1), n_total);
     launch_check();
 }
 

src/kernels/cuda/native_ple_postops.cu

@@ -118,6 +118,34 @@ __global__ void rms_rep_kernel(const float* __restrict__ input, const float* __r
     const float scale = rsqrtf(mean + NG_RMS_EPS);
     for (int col = tid; col < N; col += BlockSize) output[col] = scale * input[col] * gamma[col];
 }
+// rms_rep_kernel with the column count as a RUN-TIME divisor, as weighted_rms_norm<1024> (native_gr_rms_norm_weighted,
+// the per-token path) divides: under --use_fast_math a division by the constant N may become a multiplication by its
+// rounded reciprocal, which is not bitwise the run-time division.
+__global__ void rms_rep_rt_kernel(const float* __restrict__ input, const float* __restrict__ gamma,
+                                  float* __restrict__ output, int n_cols, float epsilon) {
+    constexpr int BlockSize = 1024;
+    const int tid = threadIdx.x;
+    const size_t row_offset = size_t(blockIdx.x) * n_cols;
+    input += row_offset;
+    output += row_offset;
+    gamma += size_t(blockIdx.x % H) * n_cols;
+    float partial = 0.0f;
+    for (int col = tid; col < n_cols; col += BlockSize) {
+        const float value = input[col];
+        partial += value * value;
+    }
+    __shared__ float sums[32];
+    partial = norm_warp_sum(partial);
+    const int lane = tid % 32;
+    if (lane == 0) sums[tid / 32] = partial;
+    __syncthreads();
+    partial = 0.0f;
+    if (lane < BlockSize / 32) partial = sums[lane];
+    partial = norm_warp_sum(partial);
+    const float mean = partial / n_cols;
+    const float scale = rsqrtf(mean + epsilon);
+    for (int col = tid; col < n_cols; col += BlockSize) output[col] = scale * input[col] * gamma[col];
+}
 __global__ void broadcast_batch_kernel(const float* value, const float* gate, float* gated, int T) {
     const size_t i = size_t(blockIdx.x) * blockDim.x + threadIdx.x;
     if (i >= size_t(T) * D) return;
@@ -156,6 +184,29 @@ __global__ void history_batch_kernel(float* history, const float* normalized, in
     for (int r = 0; r < HISTORY; ++r) history[size_t(c) * HISTORY + r] = h[r];
 }
 
+// history_batch_kernel that also writes the history after each of the T tokens (`snap`, T x 9 x 10240): pure copies,
+// the values T history advances (ple_history_advance) leave
+__global__ void history_batch_snap_kernel(float* history, const float* normalized, int T, float* snap) {
+    const int c = blockIdx.x * blockDim.x + threadIdx.x;
+    if (c >= D) return;
+    float h0[HISTORY];
+#pragma unroll
+    for (int r = 0; r < HISTORY; ++r) h0[r] = history[size_t(c) * HISTORY + r];
+    for (int t = 0; t < T; ++t) {
+        float* out = snap + size_t(t) * HISTORY * D + size_t(c) * HISTORY;
+#pragma unroll
+        for (int r = 0; r < HISTORY; ++r) {
+            const int p = r + t + 1 - HISTORY;      // after t + 1 advances row r holds token p, or an older row
+            out[r] = p >= 0 ? normalized[size_t(p) * D + c] : h0[r + t + 1];
+        }
+    }
+#pragma unroll
+    for (int r = 0; r < HISTORY; ++r) {
+        const int p = r + T - HISTORY;
+        history[size_t(c) * HISTORY + r] = p >= 0 ? normalized[size_t(p) * D + c] : h0[r + T];
+    }
+}
+
 struct Span { const void* p; size_t bytes; size_t alignment; };
 bool overlaps(Span a, Span b) {
     const auto x = reinterpret_cast<uintptr_t>(a.p), y = reinterpret_cast<uintptr_t>(b.p);
@@ -218,4 +269,21 @@ void native_ple_postops_batch(float* key, float* hidden, const float* value, flo
     history_batch_kernel<<<D / 256, 256, 0, st>>>(history, query_norm, T);
     launch_check();
 }
+
+void native_ple_postops_batch_snap(float* key, float* hidden, const float* value, float* history, const PleWeights& w,
+                                   float* query_norm, float* gated, float* gate, int T, float* snap, void* stream) {
+    if (!stream || T <= 0 || !key || !hidden || !value || !history || !query_norm || !gated || !gate || !snap)
+        throw std::invalid_argument("native PLE postops batch (snapshots): null input or empty batch");
+    auto st = static_cast<cudaStream_t>(stream);
+    const unsigned rows = unsigned(T) * H;
+    const unsigned blocks = unsigned((size_t(T) * D + 255) / 256);
+    rms_rep_rt_kernel<<<rows, 1024, 0, st>>>(key, w.norm_key, key, N, NG_RMS_EPS);
+    rms_rep_rt_kernel<<<rows, 1024, 0, st>>>(hidden, w.norm_query, query_norm, N, NG_RMS_EPS);
+    gate_kernel<<<rows, 512, 0, st>>>(key, query_norm, gate, 1.0f / std::sqrt(float(N)));
+    broadcast_batch_kernel<<<blocks, 256, 0, st>>>(value, gate, gated, T);
+    rms_rep_rt_kernel<<<rows, 1024, 0, st>>>(gated, w.norm_conv, query_norm, N, NG_RMS_EPS);
+    conv_residual_batch_kernel<<<blocks, 256, 0, st>>>(history, query_norm, w.conv1d_f16, hidden, gated, T);
+    history_batch_snap_kernel<<<D / 256, 256, 0, st>>>(history, query_norm, T, snap);
+    launch_check();
+}
 } // namespace strata::kernels

src/kernels/cuda/verify_kernels.cu

@@ -3,14 +3,21 @@
 // The per-token arithmetic of every kernel here is transcribed from its single-token original (fused_gdn.cu,
 // elementwise.cu) with the same operation order, so a verify window reproduces plain decode bit for bit.
 #include "strata/kernels/verify_kernels.hpp"
+#include "strata/core/emulate.hpp"
 #include "strata/kernels/dp4a.hpp"
+#include "strata/kernels/pdl.hpp"
 
 #include <cuda_fp16.h>
 #include <cuda_runtime.h>
 
 #include <algorithm>
+#include <atomic>
+#include <cstdint>
 #include <cstdio>
 #include <cstdlib>
+#include <mutex>
+#include <utility>
+#include <vector>
 
 namespace strata::kernels {
 namespace {
@@ -1211,8 +1218,15 @@ void copy_indexed(float* dst, const float* src, int64_t stride, const int32_t* i
     check("copy_indexed");
 }
 
-// a GPU timestamp (ns, %globaltimer) into buf[i] - the verify window's stage profiler
-namespace { __global__ void gpu_stamp_kernel(unsigned long long* buf, int i) {
+// a GPU timestamp (ns, %globaltimer) into buf[i] - the verify window's stage profiler.  Pdl: launched inside a PDL
+// stretch, the stamp lets the next kernel launch at once and takes its time when the one before it has finished.
+namespace {
+template <bool Pdl>
+__global__ void gpu_stamp_kernel(unsigned long long* buf, int i) {
+    if constexpr (Pdl) {
+        pdl_trigger();
+        pdl_wait();
+    }
     unsigned long long t;
 #if defined(STRATA_HIP_GFX906)
     t = wall_clock64() * 40ull;   // gfx906: the wall clock runs at 25 MHz (hipDeviceAttributeWallClockRate) -> ns
@@ -1222,9 +1236,121 @@ namespace { __global__ void gpu_stamp_kernel(unsigned long long* buf, int i) {
     asm volatile("mov.u64 %0, %%globaltimer;" : "=l"(t));
 #endif
     buf[i] = t;
-} }
+}
+}  // namespace
 void gpu_stamp(unsigned long long* buf, int i, void* stream) {
-    gpu_stamp_kernel<<<1, 1, 0, (cudaStream_t) stream>>>(buf, i);
+    if (pdl_scope()) launch_pdl(gpu_stamp_kernel<true>, dim3(1), dim3(1), 0, (cudaStream_t) stream, buf, i);
+    else gpu_stamp_kernel<false><<<1, 1, 0, (cudaStream_t) stream>>>(buf, i);
 }
 
+namespace {
+__global__ void copy_rows_strided_kernel(float4* __restrict__ dst, const float4* __restrict__ src, long long rows,
+                                         int w4, int src_w4) {
+    const long long total = rows * w4;
+    for (long long i = (long long) blockIdx.x * blockDim.x + threadIdx.x; i < total;
+         i += (long long) gridDim.x * blockDim.x) {
+        const long long r = i / w4, c = i - r * w4;
+        dst[i] = src[r * src_w4 + c];
+    }
+}
+}  // namespace
+void copy_rows_strided(float* dst, const float* src, int64_t rows, int64_t w, int64_t src_w, void* stream) {
+    if (rows <= 0 || w <= 0) return;
+    if ((w & 3) || (src_w & 3) || src_w < w || ((uintptr_t) dst & 15) || ((uintptr_t) src & 15)) {
+        std::fprintf(stderr, "copy_rows_strided: widths must be multiples of 4 floats (src_w >= w), pointers 16-byte "
+                             "aligned\n");
+        std::exit(1);
+    }
+    const long long total = (long long) rows * (w / 4);
+    const unsigned blocks = (unsigned) ((total + 255) / 256 < 256 ? (total + 255) / 256 : 256);
+    copy_rows_strided_kernel<<<blocks, 256, 0, (cudaStream_t) stream>>>((float4*) dst, (const float4*) src, rows,
+                                                                         (int) (w / 4), (int) (src_w / 4));
+    check("copy_rows_strided");
+}
+
+// ---- programmatic dependent launch (pdl.hpp)
+bool& pdl_scope() {
+    static thread_local bool on = false;
+    return on;
+}
+
+bool pdl_supported() {
+#if defined(__HIPCC__) || !defined(CUDART_VERSION) || CUDART_VERSION < 12030
+    return false;   // HIP; or a CUDA runtime without the capture-dependency query pdl_launch_ok needs
+#else
+    // per device (a layer split runs on several): 1 = runs, 2 = does not.  sm_90+ (the card's, or STRATA_EMULATE_CC's),
+    // and code built for it: on a build with only older code the driver JIT-compiles PTX without griddepcontrol.
+    static std::atomic<int> state[64];
+    int dev = 0;
+    if (cudaGetDevice(&dev) != cudaSuccess || dev < 0 || dev >= 64) { cudaGetLastError(); return false; }
+    int s = state[dev].load(std::memory_order_relaxed);
+    if (s == 0) {
+        static const bool env_on = [] { const char* v = std::getenv("STRATA_DF_PDL"); return v != nullptr && std::atoi(v) != 0; }();   // opt-in until the A/B says otherwise
+        int major = 0;
+        cudaFuncAttributes fa{};
+        const bool on = env_on && cudaDeviceGetAttribute(&major, cudaDevAttrComputeCapabilityMajor, dev) == cudaSuccess &&
+                        strata::cc_major_of(major) >= 9 &&
+                        cudaFuncGetAttributes(&fa, gpu_stamp_kernel<true>) == cudaSuccess && fa.ptxVersion >= 90 &&
+                        fa.binaryVersion >= 90;
+        cudaGetLastError();
+        s = on ? 1 : 2;
+        state[dev].store(s, std::memory_order_relaxed);
+    }
+    return s == 1;
+#endif
+}
+
+#if !defined(__HIPCC__)
+bool pdl_launch_ok(const void* kernel, cudaStream_t stream) {
+#if !defined(CUDART_VERSION) || CUDART_VERSION < 12030
+    (void) kernel; (void) stream;
+    return false;
+#else
+    // the kernel's own code is sm_90+ (each kernel checked once per device: its PTX, or a JIT from older PTX, decides)
+    {
+        static std::mutex mu;
+        static std::vector<std::pair<std::pair<const void*, int>, bool>> seen;
+        int dev = 0;
+        if (cudaGetDevice(&dev) != cudaSuccess) { cudaGetLastError(); return false; }
+        std::lock_guard<std::mutex> lock(mu);
+        bool found = false, ok = false;
+        for (const auto& e : seen)
+            if (e.first.first == kernel && e.first.second == dev) { found = true; ok = e.second; break; }
+        if (!found) {
+            cudaFuncAttributes fa{};
+            ok = cudaFuncGetAttributes(&fa, kernel) == cudaSuccess && fa.ptxVersion >= 90 && fa.binaryVersion >= 90;
+            cudaGetLastError();
+            seen.push_back({{kernel, dev}, ok});
+        }
+        if (!ok) return false;
+    }
+    // every node the next captured node will depend on is a kernel (CUDA allows a programmatic edge only between two
+    // kernel nodes): never after a memcpy, memset, host or event node.  A join of branches gives several kernels, a
+    // programmatic edge from each; STRATA_DF_PDL=2 allows a single predecessor only.
+    static const int mode = [] { const char* v = std::getenv("STRATA_DF_PDL"); return v ? std::atoi(v) : 1; }();
+    cudaStreamCaptureStatus status = cudaStreamCaptureStatusNone;
+    const cudaGraphNode_t* deps = nullptr;
+    const cudaGraphEdgeData* edges = nullptr;
+    size_t n = 0;
+#if CUDART_VERSION >= 13000
+    const cudaError_t e = cudaStreamGetCaptureInfo(stream, &status, nullptr, nullptr, &deps, &edges, &n);
+#else
+    const cudaError_t e = cudaStreamGetCaptureInfo_v3(stream, &status, nullptr, nullptr, &deps, &edges, &n);
+#endif
+    if (e != cudaSuccess) { cudaGetLastError(); return false; }
+    if (status != cudaStreamCaptureStatusActive || n == 0 || deps == nullptr) return false;
+    if (mode == 2 && n != 1) return false;
+    for (size_t i = 0; i < n; ++i) {
+        if (edges != nullptr && (edges[i].type != cudaGraphDependencyTypeDefault || edges[i].from_port != 0 ||
+                                 edges[i].to_port != 0))
+            return false;
+        cudaGraphNodeType type;
+        if (cudaGraphNodeGetType(deps[i], &type) != cudaSuccess) { cudaGetLastError(); return false; }
+        if (type != cudaGraphNodeTypeKernel) return false;
+    }
+    return true;
+#endif
+}
+#endif
+
 }  // namespace strata::kernels

src/kernels/pdl_parity.cpp

@@ -0,0 +1,306 @@
+// src/kernels/pdl_parity.cpp - programmatic dependent launch (pdl.hpp) in a captured graph, tested bitwise.
+//
+//     build/pdl_parity [--replays N]
+//
+// WHAT IT CHECKS.  A chain shaped like the verify window's - quantize, multi-column MMVQ (a small-K and a large-K
+// format), BF16 multi-row GEMV, the shared Q8_1 buffer rewritten by the next quantization, a memcpy node, a fork onto a
+// second stream and a join - is captured twice: with pdl_scope() off (ordinary launches) and on.  Then:
+//   1. every output of the PDL graph is bitwise the ordinary graph's, over many replays; before each replay every
+//      intermediate is overwritten with NaN bytes, so a kernel that read its input before the kernel writing it had
+//      finished would see NaNs and fail;
+//   2. every programmatic edge of the PDL graph joins two kernel nodes, and every node with a programmatic in-edge
+//      has kernels only before it (the kernel after the memcpy must be an ordinary launch); with STRATA_DF_PDL=2 also
+//      no programmatic edge into a node with a second in-edge (the kernel after the join);
+//   3. power: where pdl_supported() (sm_90+ with code built for it) the PDL graph has programmatic edges, elsewhere
+//      it has none.
+// Also copy_rows_strided against the cudaMemcpy2DAsync it replaces in the verify window.  GPU, synthetic, no model.
+#include "strata/kernels/bf16_bits.hpp"
+#include "strata/kernels/bf16_gemv.hpp"
+#include "strata/kernels/native_mmvq.hpp"
+#include "strata/kernels/pdl.hpp"
+#include "strata/kernels/verify_kernels.hpp"
+
+#include <cuda_runtime.h>
+
+#include <cmath>
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <map>
+#include <random>
+#include <string>
+#include <vector>
+
+namespace k = strata::kernels;
+
+namespace {
+
+int g_fail = 0;
+
+void ck(cudaError_t e, const char* what) {
+    if (e != cudaSuccess) {
+        std::printf("CUDA: %s: %s\n", what, cudaGetErrorString(e));
+        std::exit(2);
+    }
+}
+
+// The set-up fills and copies run on the legacy default stream, which does not order the non-blocking streams the
+// kernels run on, and may return before the device has finished them: each one is completed before it is used.
+template <typename T>
+T* dalloc(size_t n) {
+    T* p = nullptr;
+    ck(cudaMalloc(&p, n * sizeof(T)), "malloc");
+    ck(cudaMemset(p, 0, n * sizeof(T)), "memset");
+    ck(cudaDeviceSynchronize(), "memset sync");
+    return p;
+}
+
+// a normal fp16 in +-[2^-10, 2^-5) (mmvq_multi_parity's: finite scales keep every output finite)
+uint16_t sane_half(std::mt19937& rng) {
+    const uint32_t r = rng();
+    return (uint16_t) (((r >> 31) << 15) | ((5u + (r >> 10) % 5u) << 10) | (r & 0x3ffu));
+}
+
+// random quantized weights with sane fp16 scales at byte 0 of each block (Q8_0: 34 bytes, IQ4_XS: 136 bytes)
+void* quant_weights(int type, int n_in, int n_out, int block_elems, int block_bytes, std::mt19937& rng) {
+    const size_t bytes = k::native_mmvq_weight_bytes(type, n_in, n_out);
+    std::vector<uint8_t> w(bytes);
+    for (auto& b : w) b = (uint8_t) (rng() & 0xff);
+    const size_t blocks = (size_t) n_out * (size_t) (n_in / block_elems);
+    if (blocks * (size_t) block_bytes != bytes) {
+        std::printf("type %d: %zu blocks of %d bytes are not the %zu weight bytes\n", type, blocks, block_bytes, bytes);
+        std::exit(2);
+    }
+    for (size_t i = 0; i < blocks; ++i) {
+        const uint16_t h = sane_half(rng);
+        std::memcpy(&w[i * (size_t) block_bytes], &h, 2);
+    }
+    void* d = nullptr;
+    ck(cudaMalloc(&d, bytes), "malloc w");
+    ck(cudaMemcpy(d, w.data(), bytes, cudaMemcpyHostToDevice), "copy w");
+    ck(cudaDeviceSynchronize(), "copy w sync");
+    return d;
+}
+
+uint16_t* bf16_weights(int n_in, int n_out, float scale, std::mt19937& rng) {
+    std::normal_distribution<float> nd(0.f, scale);
+    std::vector<uint16_t> w((size_t) n_in * n_out);
+    for (auto& v : w) v = k::bf16_from_f32(nd(rng));
+    uint16_t* d = dalloc<uint16_t>(w.size());
+    ck(cudaMemcpy(d, w.data(), w.size() * 2, cudaMemcpyHostToDevice), "copy bf16 w");
+    ck(cudaDeviceSynchronize(), "copy bf16 w sync");
+    return d;
+}
+
+constexpr int T = 4, N = 2560, NB = 6144;   // a window of 4 tokens, n_embd, a wide projection
+
+struct Bufs {
+    float *x0, *y1, *y2, *y3, *y4, *y5, *y5c, *y6, *y7, *y8;
+    uint8_t *q, *q2;
+};
+
+// The chain (one capture).  `side`/`ev*`: the fork and join, as the verify window's branches make them.
+void record(const Bufs& b, void* w_q8, void* w_iq4, uint16_t* w_b1, uint16_t* w_b2, cudaStream_t s, cudaStream_t side,
+            cudaEvent_t ev_fork, cudaEvent_t ev_join) {
+    k::native_quantize_q8_1(b.x0, b.q, N, T, s);                          // A
+    k::native_mmvq(8, w_q8, b.q, b.y1, N, N, T, s);                       // B  Q8_0, large K
+    k::native_quantize_q8_1(b.y1, b.q, N, T, s);                          // C  rewrites the q B read
+    k::native_mmvq(23, w_iq4, b.q, b.y2, N, NB, T, s);                    // D  IQ4_XS, small K
+    k::bf16_gemv_fp32_mmvf_multi(b.y2, NB, w_b1, b.y3, N, NB, N, T, s);  // E
+    ck(cudaMemcpyAsync(b.y4, b.y3, (size_t) T * N * 4, cudaMemcpyDeviceToDevice, s), "memcpy node");   // F
+    k::native_quantize_q8_1(b.y4, b.q, N, T, s);                          // G  after a memcpy: ordinary
+    k::native_mmvq(8, w_q8, b.q, b.y5, N, N, T, s);                       // H
+    ck(cudaEventRecord(ev_fork, s), "fork");
+    ck(cudaStreamWaitEvent(side, ev_fork, 0), "fork wait");
+    k::copy_rows_strided(b.y5c, b.y5, T, N, N, side);                     // side: an ordinary kernel after H
+    k::native_quantize_q8_1(b.y5c, b.q2, N, T, side);
+    k::native_mmvq(8, w_q8, b.q2, b.y6, N, N, T, side);
+    k::bf16_gemv_fp32_mmvf_multi(b.y5, N, w_b2, b.y7, N, N, N, T, s);    // main: beside the side branch
+    ck(cudaEventRecord(ev_join, side), "join");
+    ck(cudaStreamWaitEvent(s, ev_join, 0), "join wait");
+    k::bf16_gemv_fp32_mmvf_multi(b.y6, N, w_b2, b.y8, N, N, N, T, s);    // after the join: ordinary
+}
+
+cudaGraphExec_t capture(bool pdl, const Bufs& b, void* w_q8, void* w_iq4, uint16_t* w_b1, uint16_t* w_b2,
+                        cudaStream_t s, cudaStream_t side, cudaEvent_t ef, cudaEvent_t ej, size_t& prog, size_t& bad,
+                        size_t& multi, size_t& mixed) {
+    k::pdl_scope() = pdl;
+    ck(cudaStreamBeginCapture(s, cudaStreamCaptureModeThreadLocal), "begin capture");
+    record(b, w_q8, w_iq4, w_b1, w_b2, s, side, ef, ej);
+    cudaGraph_t g = nullptr;
+    ck(cudaStreamEndCapture(s, &g), "end capture");
+    k::pdl_scope() = false;
+    prog = bad = multi = mixed = 0;
+#if CUDART_VERSION >= 12030
+    size_t ne = 0;
+#if CUDART_VERSION >= 13000
+    ck(cudaGraphGetEdges(g, nullptr, nullptr, nullptr, &ne), "edges");
+#else
+    ck(cudaGraphGetEdges_v2(g, nullptr, nullptr, nullptr, &ne), "edges");
+#endif
+    std::vector<cudaGraphNode_t> from(ne), to(ne);
+    std::vector<cudaGraphEdgeData> ed(ne);
+#if CUDART_VERSION >= 13000
+    ck(cudaGraphGetEdges(g, from.data(), to.data(), ed.data(), &ne), "edges");
+#else
+    ck(cudaGraphGetEdges_v2(g, from.data(), to.data(), ed.data(), &ne), "edges");
+#endif
+    std::map<cudaGraphNode_t, int> in_edges, prog_in, nonkernel_in;
+    for (size_t i = 0; i < ne; ++i) {
+        ++in_edges[to[i]];
+        cudaGraphNodeType ta;
+        ck(cudaGraphNodeGetType(from[i], &ta), "node type");
+        if (ta != cudaGraphNodeTypeKernel) ++nonkernel_in[to[i]];
+        if (ed[i].type == cudaGraphDependencyTypeProgrammatic) ++prog_in[to[i]];
+    }
+    for (size_t i = 0; i < ne; ++i) {
+        if (ed[i].type != cudaGraphDependencyTypeProgrammatic) continue;
+        ++prog;
+        cudaGraphNodeType ta, tb;
+        ck(cudaGraphNodeGetType(from[i], &ta), "node type");
+        ck(cudaGraphNodeGetType(to[i], &tb), "node type");
+        if (ta != cudaGraphNodeTypeKernel || tb != cudaGraphNodeTypeKernel) ++bad;
+        if (in_edges[to[i]] != 1) ++multi;
+    }
+    for (const auto& kv : prog_in)
+        if (nonkernel_in[kv.first] != 0) ++mixed;
+#endif
+    cudaGraphExec_t ex = nullptr;
+    ck(cudaGraphInstantiate(&ex, g, 0), "instantiate");
+    ck(cudaGraphDestroy(g), "destroy graph");
+    return ex;
+}
+
+std::vector<uint32_t> read(const float* d, size_t n) {
+    std::vector<uint32_t> h(n);
+    ck(cudaMemcpy(h.data(), d, n * 4, cudaMemcpyDeviceToHost), "read");
+    return h;
+}
+
+void copy_rows_strided_check(cudaStream_t s) {
+    const int64_t rows = 3 * 32, w = 256, src_w = 512;   // a 3-token window's q heads: query half of each q/gate pair
+    std::mt19937 rng(7);
+    std::normal_distribution<float> nd(0.f, 1.f);
+    std::vector<float> src((size_t) rows * src_w);
+    for (auto& v : src) v = nd(rng);
+    float* ds = dalloc<float>(src.size());
+    float* da = dalloc<float>((size_t) rows * w);
+    float* db = dalloc<float>((size_t) rows * w);
+    ck(cudaMemcpy(ds, src.data(), src.size() * 4, cudaMemcpyHostToDevice), "copy src");
+    ck(cudaDeviceSynchronize(), "copy src sync");
+    k::copy_rows_strided(da, ds, rows, w, src_w, s);
+    ck(cudaMemcpy2DAsync(db, (size_t) w * 4, ds, (size_t) src_w * 4, (size_t) w * 4, (size_t) rows,
+                         cudaMemcpyDeviceToDevice, s), "memcpy2d");
+    ck(cudaStreamSynchronize(s), "sync");
+    const auto a = read(da, (size_t) rows * w), b = read(db, (size_t) rows * w);
+    const bool ok = a == b;
+    std::printf("copy_rows_strided vs cudaMemcpy2DAsync (%lld rows of %lld of %lld): %s\n", (long long) rows,
+                (long long) w, (long long) src_w, ok ? "bitwise equal" : "DIFFERENT");
+    if (!ok) ++g_fail;
+    cudaFree(ds); cudaFree(da); cudaFree(db);
+}
+
+}  // namespace
+
+int main(int argc, char** argv) {
+    int replays = 200;
+    for (int i = 1; i < argc; ++i)
+        if (!std::strcmp(argv[i], "--replays") && i + 1 < argc) replays = std::atoi(argv[++i]);
+    int dev = 0;
+    cudaDeviceProp prop{};
+    ck(cudaGetDevice(&dev), "device");
+    ck(cudaGetDeviceProperties(&prop, dev), "props");
+    const bool supported = k::pdl_supported();
+    std::printf("pdl_parity: %s (sm_%d%d), pdl_supported() = %d\n", prop.name, prop.major, prop.minor, (int) supported);
+
+    cudaStream_t s = nullptr, side = nullptr;
+    cudaEvent_t ef = nullptr, ej = nullptr;
+    ck(cudaStreamCreateWithFlags(&s, cudaStreamNonBlocking), "stream");
+    ck(cudaStreamCreateWithFlags(&side, cudaStreamNonBlocking), "stream");
+    ck(cudaEventCreateWithFlags(&ef, cudaEventDisableTiming), "event");
+    ck(cudaEventCreateWithFlags(&ej, cudaEventDisableTiming), "event");
+
+    copy_rows_strided_check(s);
+
+    std::mt19937 rng(1234);
+    void* w_q8 = quant_weights(8, N, N, 32, 34, rng);
+    void* w_iq4 = quant_weights(23, N, NB, 256, 136, rng);
+    uint16_t* w_b1 = bf16_weights(NB, N, 0.02f, rng);
+    uint16_t* w_b2 = bf16_weights(N, N, 0.02f, rng);
+    Bufs b{};
+    b.x0 = dalloc<float>((size_t) T * N);
+    float** outs[] = {&b.y1, &b.y2, &b.y3, &b.y4, &b.y5, &b.y5c, &b.y6, &b.y7, &b.y8};
+    constexpr int NOUT = 9;
+    const size_t out_n[NOUT] = {(size_t) T * N, (size_t) T * NB, (size_t) T * N, (size_t) T * N, (size_t) T * N,
+                                (size_t) T * N, (size_t) T * N, (size_t) T * N, (size_t) T * N};
+    for (int i = 0; i < NOUT; ++i) *outs[i] = dalloc<float>(out_n[i]);
+    const size_t qbytes = k::native_q8_1_bytes(N, T);
+    b.q = dalloc<uint8_t>(qbytes);
+    b.q2 = dalloc<uint8_t>(qbytes);
+    {
+        std::normal_distribution<float> nd(0.f, 1.f);
+        std::vector<float> x((size_t) T * N);
+        for (auto& v : x) v = nd(rng);
+        ck(cudaMemcpy(b.x0, x.data(), x.size() * 4, cudaMemcpyHostToDevice), "copy x");
+        ck(cudaDeviceSynchronize(), "copy x sync");
+    }
+
+    size_t prog0, bad0, multi0, mixed0, prog1, bad1, multi1, mixed1;
+    cudaGraphExec_t plain = capture(false, b, w_q8, w_iq4, w_b1, w_b2, s, side, ef, ej, prog0, bad0, multi0, mixed0);
+    cudaGraphExec_t pdl = capture(true, b, w_q8, w_iq4, w_b1, w_b2, s, side, ef, ej, prog1, bad1, multi1, mixed1);
+    const char* mode = std::getenv("STRATA_DF_PDL");
+    const bool single = mode != nullptr && std::atoi(mode) == 2;
+    std::printf("graph edges: ordinary capture %zu programmatic; PDL capture %zu programmatic, %zu not kernel to kernel, "
+                "%zu nodes with a programmatic and a non-kernel in-edge, %zu into a node with other in-edges%s\n", prog0,
+                prog1, bad1, mixed1, multi1, single ? " (STRATA_DF_PDL=2: must be 0)" : "");
+    if (prog0 != 0 || bad1 != 0 || mixed1 != 0 || (single && multi1 != 0)) ++g_fail;
+#if CUDART_VERSION >= 12030
+    if (supported && prog1 == 0) {
+        std::printf("FAIL: PDL is supported here but the PDL capture has no programmatic edge (the test has no power)\n");
+        ++g_fail;
+    }
+    if (!supported && prog1 != 0) {
+        std::printf("FAIL: PDL is not supported here but the capture has programmatic edges\n");
+        ++g_fail;
+    }
+#endif
+
+    auto poison = [&]() {   // NaN bytes in every intermediate and output
+        for (int i = 0; i < NOUT; ++i) ck(cudaMemsetAsync(*outs[i], 0xff, out_n[i] * 4, s), "poison");
+        ck(cudaMemsetAsync(b.q, 0xff, qbytes, s), "poison");
+        ck(cudaMemsetAsync(b.q2, 0xff, qbytes, s), "poison");
+    };
+    poison();
+    ck(cudaGraphLaunch(plain, s), "launch plain");
+    ck(cudaStreamSynchronize(s), "sync plain");
+    std::vector<std::vector<uint32_t>> ref;
+    for (int i = 0; i < NOUT; ++i) ref.push_back(read(*outs[i], out_n[i]));
+    long long nonfinite = 0;
+    for (const auto& r : ref)
+        for (uint32_t v : r) {
+            float f;
+            std::memcpy(&f, &v, 4);
+            if (!std::isfinite(f)) ++nonfinite;
+        }
+    if (nonfinite) {
+        std::printf("FAIL: the ordinary graph left %lld non-finite outputs (the data has no power)\n", nonfinite);
+        ++g_fail;
+    }
+    int bad_replays = 0;
+    for (int r = 0; r < replays; ++r) {
+        poison();
+        ck(cudaGraphLaunch(pdl, s), "launch pdl");
+        ck(cudaStreamSynchronize(s), "sync pdl");
+        bool same = true;
+        for (int i = 0; i < NOUT && same; ++i) same = read(*outs[i], out_n[i]) == ref[(size_t) i];
+        if (!same) ++bad_replays;
+    }
+    std::printf("PDL graph vs ordinary graph over %d replays: %d differ\n", replays, bad_replays);
+    if (bad_replays) ++g_fail;
+
+    cudaGraphExecDestroy(plain);
+    cudaGraphExecDestroy(pdl);
+    std::printf("pdl_parity: %s\n", g_fail ? "FAIL" : "PASS");
+    return g_fail ? 1 : 0;
+}

src/kernels/verify_batch_parity.cpp

@@ -0,0 +1,173 @@
+// src/kernels/verify_batch_parity.cpp - the verify window's one-launch forms (the hyper-connection write, the PLE post-ops) against the per-token launches they
+// replace, bitwise.
+//
+//     build/verify_batch_parity
+//
+// 1. native_gr_post_multi against n native_gr_post calls (the hyper-connection write before the PLE block).
+// 2. native_ple_postops_batch_snap (STRATA_DF_PLE) against the per-token loop of the verify window: per token
+//    native_ple_postops, ple_history_advance and a copy of the history (the commit's snapshot); the results, the final
+//    history and every snapshot, for windows of 2..8 tokens.
+// GPU, synthetic, no model.
+#include "strata/kernels/f16_bits.hpp"
+#include "strata/kernels/kv_q4.hpp"
+#include "strata/kernels/kv_stream.hpp"
+#include "strata/kernels/native_gr_postops.hpp"
+#include "strata/kernels/native_ple_postops.hpp"
+#include "strata/kernels/ngram.hpp"
+#include "strata/kernels/ple.hpp"
+#include "strata/kernels/qsa.hpp"
+
+#include <cuda_runtime.h>
+
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <random>
+#include <vector>
+
+namespace k = strata::kernels;
+
+namespace {
+
+int g_fail = 0;
+
+void ck(cudaError_t e, const char* what) {
+    if (e != cudaSuccess) {
+        std::printf("CUDA: %s: %s\n", what, cudaGetErrorString(e));
+        std::exit(2);
+    }
+}
+
+// The set-up copies and fills run on the legacy default stream, which does NOT order work on the non-blocking stream
+// the kernels run on (and cudaMemset / a pageable cudaMemcpy may return before the device has finished them): each
+// one is completed here, or a kernel could read its inputs - or have its outputs overwritten - while they land.
+template <typename T>
+T* dalloc(size_t n, int fill = 0) {
+    T* p = nullptr;
+    ck(cudaMalloc(&p, n * sizeof(T) + 64), "malloc");
+    ck(cudaMemset(p, fill, n * sizeof(T) + 64), "memset");
+    ck(cudaDeviceSynchronize(), "memset sync");
+    return p;
+}
+
+template <typename T>
+std::vector<T> read(const T* d, size_t n) {
+    std::vector<T> h(n);
+    ck(cudaMemcpy(h.data(), d, n * sizeof(T), cudaMemcpyDeviceToHost), "read");
+    return h;
+}
+
+template <typename T>
+T* upload(const std::vector<T>& h) {
+    T* d = dalloc<T>(h.size());
+    ck(cudaMemcpy(d, h.data(), h.size() * sizeof(T), cudaMemcpyHostToDevice), "upload");
+    ck(cudaDeviceSynchronize(), "upload sync");
+    return d;
+}
+
+void report(const char* what, bool ok) {
+    std::printf("%-64s %s\n", what, ok ? "bitwise equal" : "DIFFERENT");
+    if (!ok) ++g_fail;
+}
+
+// ---- 1. the hyper-connection write
+void gr_post(cudaStream_t s) {
+    const int N = 2560, HC = 4;
+    std::mt19937 rng(12);
+    std::normal_distribution<float> nd(0.f, 1.f);
+    for (int n : {2, 4, 8}) {
+        std::vector<float> r((size_t) n * HC * N), bo((size_t) n * N), inj((size_t) n * HC);
+        for (auto& x : r) x = nd(rng);
+        for (auto& x : bo) x = nd(rng);
+        for (auto& x : inj) x = 3.f * nd(rng);
+        float* r0 = upload(r);
+        float* r1 = upload(r);
+        float* d_bo = upload(bo);
+        float* d_inj = upload(inj);
+        for (int t = 0; t < n; ++t)
+            k::native_gr_post(r0 + (size_t) t * HC * N, d_bo + (size_t) t * N, d_inj + (size_t) t * HC,
+                              r0 + (size_t) t * HC * N, N, HC, s);
+        k::native_gr_post_multi(r1, d_bo, d_inj, r1, N, HC, n, (long long) HC * N, N, HC, s);
+        ck(cudaStreamSynchronize(s), "gr post");
+        char what[128];
+        std::snprintf(what, sizeof what, "native_gr_post_multi, %d tokens (in place)", n);
+        report(what, read(r0, r.size()) == read(r1, r.size()));
+        cudaFree(r0); cudaFree(r1); cudaFree(d_bo); cudaFree(d_inj);
+    }
+}
+
+// ---- 2. the PLE post-ops of a window, with the commit's history snapshots
+void ple(cudaStream_t s) {
+    constexpr int N = k::NG_N_EMBD, D = k::NG_HC_DIM, HS = k::NG_HIST * k::NG_HC_DIM;
+    std::mt19937 rng(13);
+    std::normal_distribution<float> nd(0.f, 1.f);
+    std::vector<float> nk(D), nq(D), nc(D);
+    for (auto* v : {&nk, &nq, &nc})
+        for (auto& x : *v) x = 1.f + 0.1f * nd(rng);
+    std::vector<uint16_t> conv((size_t) 4 * D);
+    for (auto& x : conv) x = k::f16_from_f32(0.3f * nd(rng));
+    k::PleWeights w;
+    w.norm_key = upload(nk);
+    w.norm_query = upload(nq);
+    w.norm_conv = upload(nc);
+    w.conv1d_f16 = upload(conv);
+    std::vector<float> hist0((size_t) HS);
+    for (auto& x : hist0) x = nd(rng);
+    for (int n : {2, 3, 4, 8}) {
+        std::vector<float> key((size_t) n * D), hid((size_t) n * D), val((size_t) n * N);
+        for (auto& x : key) x = 4.f * nd(rng);
+        for (auto& x : hid) x = nd(rng);
+        for (auto& x : val) x = nd(rng);
+        // per token, as the verify window's loop: native_ple_postops (ple_block_projected's buffers), the history
+        // advance and the snapshot copy
+        float* key0 = upload(key);
+        float* hid0 = upload(hid);
+        float* val0 = upload(val);
+        float* hist_a = upload(hist0);
+        float* snap_a = dalloc<float>((size_t) n * HS);
+        float *q = dalloc<float>(D), *nrm = dalloc<float>(D), *g = dalloc<float>(4), *gd = dalloc<float>(D);
+        float* cv = dalloc<float>(D);
+        for (int t = 0; t < n; ++t) {
+            k::NativePlePostopsBuffers b{q, nrm, g, gd, nrm, cv, hid0 + (size_t) t * D};
+            k::native_ple_postops(key0 + (size_t) t * D, hid0 + (size_t) t * D, val0 + (size_t) t * N, hist_a, w, b, s);
+            k::ple_history_advance(hist_a, nrm, s);
+            ck(cudaMemcpyAsync(snap_a + (size_t) t * HS, hist_a, (size_t) HS * 4, cudaMemcpyDeviceToDevice, s), "snap");
+        }
+        // the window's one pass
+        float* key1 = upload(key);
+        float* hid1 = upload(hid);
+        float* val1 = upload(val);
+        float* hist_b = upload(hist0);
+        float* snap_b = dalloc<float>((size_t) n * HS);
+        float *qn = dalloc<float>((size_t) n * D), *gated = dalloc<float>((size_t) n * D), *gate = dalloc<float>((size_t) n * 4);
+        k::native_ple_postops_batch_snap(key1, hid1, val1, hist_b, w, qn, gated, gate, n, snap_b, s);
+        ck(cudaStreamSynchronize(s), "ple");
+        char what[128];
+        std::snprintf(what, sizeof what, "native_ple_postops_batch_snap, %d tokens: results", n);
+        report(what, read(hid0, hid.size()) == read(hid1, hid.size()));
+        std::snprintf(what, sizeof what, "native_ple_postops_batch_snap, %d tokens: history", n);
+        report(what, read(hist_a, (size_t) HS) == read(hist_b, (size_t) HS));
+        std::snprintf(what, sizeof what, "native_ple_postops_batch_snap, %d tokens: snapshots", n);
+        report(what, read(snap_a, (size_t) n * HS) == read(snap_b, (size_t) n * HS));
+        for (float* p : {key0, hid0, val0, hist_a, snap_a, q, nrm, g, gd, cv, key1, hid1, val1, hist_b, snap_b, qn, gated,
+                         gate})
+            cudaFree(p);
+    }
+}
+
+}  // namespace
+
+int main() {
+    cudaDeviceProp prop{};
+    int dev = 0;
+    ck(cudaGetDevice(&dev), "device");
+    ck(cudaGetDeviceProperties(&prop, dev), "props");
+    std::printf("verify_batch_parity: %s (sm_%d%d)\n", prop.name, prop.major, prop.minor);
+    cudaStream_t s = nullptr;
+    ck(cudaStreamCreateWithFlags(&s, cudaStreamNonBlocking), "stream");
+    gr_post(s);
+    ple(s);
+    std::printf("verify_batch_parity: %s\n", g_fail ? "FAIL" : "PASS");
+    return g_fail ? 1 : 0;
+}