Verify window: the projections' rows a warp by format and tokens, the GDN qkv grid first, IQ3_S decoded once; bitwise

- native_mmvq_il picks the rows a warp (1, 2 or 4) from a table by format, token count (1-4) and row class, measured
  alone at the engine's clocks (scratch/mm1_bench.cu) and in the engine, where the kernels run beside other branches.
  It took 2 rows, 4 for Q4_K and Q6_K from 8192 rows: at 3-4 tokens those 4-row instances need 115-193 registers, and
  Q4_K's qkv took 34-39 us where 2 rows take 25-30; Q6_K runs fastest at 2-3 tokens with 1 row (its out projection at
  3 tokens 34.6 -> 20.8 us); UD's Q8_0 head at 2 tokens with 1 row (837 -> 773 us).
- One column also goes through native_mmvq_il (its kernels read the q8_1 blocks, Q81One): IQ3_S takes them at every
  token count, and the shared expert's IQ3_S gate and up projections at 2+ tokens decode their weights once for all
  tokens (shared_expert_multi quantizes with the interleaved copy, the verify window's xil_): 11.1 / 15.1 / 17.9 ->
  9.0 / 10.5 / 8.2 us at 2 / 3 / 4 tokens.  IQ3SCols loads a block's 8 activation positions before its sums: its
  3-column instance, allotted 40 registers, loaded each after the last one's use.
- The GDN mixer records qkv's grid before the gate projection's: the conv runs beside the gate's last blocks.
- The interleaved kernels keep a grid-stride loop though the grid gives each warp one group: without it ptxas gives
  several instances fewer registers, and Q6_K's at 3 tokens run 10-50% slower.
- mmvq_il_parity: 1-8 columns with each rows a warp; --bench at the engine's clocks (a one-thread spin before each
  launch: a continuous load reaches the 3090's 260 W cap).

The mixers (IQ3_XXS thinking, us at 1 / 2 / 3 / 4 tokens): GDN 64.0 / 78.0 / 90.2 / 87.3 -> 62.7 / 71.9 / 76.9 /
84.1, QSA 105.4 / 118.9 / 144.3 / 141.3 -> 106.0 / 116.2 / 133.4 / 143.1; the shared expert ends last before the
combine in 1.3 / 0.1 / 0.8 / 1.1% of the layers (1.3 / 0.4 / 16.5 / 14.6).  Window hashes identical (thinking, code
and after 16K; both models; --spec-split and --second-gpu -1).  Fixed text, both GPUs, against aa209f7, ms per round
(6 runs of each build, half of them first in their pair: a run's place moved code-text runs by 3-6%):
IQ3_XXS code -0.9%, thinking -1.0%, after 16K +0.3%; UD-Q4_K_XL code -0.9%, thinking -1.0%, after 16K -1.3%.

Claude-Session: https://claude.ai/code/session_019XxRxQQNyk7Ycdr4ComSs1

include/strata/kernels/iq_kernels.hpp

@@ -22,9 +22,11 @@ void quantize_q8_1_rows(const float* x, int64_t n_rows, int64_t n_cols, void* y,
 
 /// y[c][r] = W[r] . x[c] for `ncols` columns of q8_1 activations (x stride n_in/32 blocks per column).
 void iq_mmvq(int ggml_type, const void* w, const void* x_q8_1, float* y, int n_in, int n_out, int ncols, void* stream);
-/// iq_mmvq for 2..8 columns from their interleaved copy (native_quantize_q8_1_il), bitwise: IQ3_S.
+/// iq_mmvq for 2..8 columns from their interleaved copy (native_quantize_q8_1_il), or one from its q8_1 blocks,
+/// bitwise, a warp taking `rows` rows (1, 2 or 4): IQ3_S.
 bool iq_mmvq_il_supported(int ggml_type) noexcept;
-void iq_mmvq_il(int ggml_type, const void* w, const void* x_il, float* y, int n_in, int n_out, int ncols, void* stream);
+void iq_mmvq_il(int ggml_type, const void* w, const void* x_q8_1, const void* x_il, float* y, int n_in, int n_out,
+                int ncols, void* stream, int rows);
 
 /// Dequantize `n` contiguous values (n a multiple of 256) to fp16 / fp32.
 void iq_dequant_f16(int ggml_type, const void* src, int64_t n, uint16_t* dst, void* stream);

include/strata/kernels/native_mmvq.hpp

@@ -116,10 +116,13 @@ constexpr int native_q8_1_il_cp(int ncols) { return ncols <= 2 ? 2 : ncols <= 4
 std::size_t native_q8_1_il_bytes(int n_in, int ncols);
 // native_quantize_q8_1 that also writes the interleaved copy to x_il (native_q8_1_il_bytes); 2 <= ncols <= 8.
 void native_quantize_q8_1_il(const float* x, void* x_q8_1, void* x_il, int n_in, int ncols, void* stream);
-// native_mmvq for 2..8 columns quantized by native_quantize_q8_1_il, bitwise native_mmvq's output: kernels that read
-// the interleaved copy for IQ4_XS, Q4_K, Q5_K, Q6_K, Q8_0 and IQ3_S matrices of at least 2048 rows, native_mmvq
-// otherwise.
+// native_mmvq for 1 column (x_il unused) or 2..8 quantized by native_quantize_q8_1_il, bitwise native_mmvq's output:
+// kernels in which a warp takes 1, 2 or 4 rows, reading 2+ columns from the interleaved copy, where they are the faster
+// on the RTX 3090 (a table by format, columns and rows: IQ4_XS, Q4_K, Q5_K, Q6_K and Q8_0 matrices of 2+ columns and
+// 2048+ rows, IQ3_S ones); native_mmvq otherwise.
 void native_mmvq_il(int ggml_type, const void* weights, const void* x_q8_1, const void* x_il, float* y, int n_in,
                     int n_out, int ncols, void* stream);
+// Tests and benchmarks: every native_mmvq_il call takes `rows` a warp (0: the table), and with `one` 1 column too.
+void native_mmvq_il_tune(int rows, bool one);
 
 } // namespace strata::kernels

include/strata/kernels/shared_expert.hpp

@@ -38,6 +38,7 @@ struct NativeSharedWeights {
     const void* up_data = nullptr;
     const void* down_data = nullptr;
     void* q8_1 = nullptr;
+    void* q8_1_il = nullptr;   // shared_expert_multi: the interleaved copy of 2+ tokens' inputs (native_q8_1_il_bytes)
 };
 
 /// Bytes of caller-owned scratch `shared_expert` needs.  **THE KERNEL USED TO `cudaMalloc` FOUR BUFFERS ON
@@ -78,7 +79,8 @@ void shared_expert_gate_rows(float* out, float* g, int64_t n_embd, int n_tok, vo
 /// scratch; `out` (n_tok, n_embd) receives the down projection's rows and `g` the scalar gate's n_tok logits,
 /// computed on `gate_stream` (may be `stream`), so that the chain ends with the down projection:
 /// `shared_expert_gate_rows` on them, or `native_moe_gather_combine` applying them, gives every token bitwise
-/// `shared_expert` on that token.  `nw.q8_1` must hold n_tok columns of n_embd.
+/// `shared_expert` on that token.  `nw.q8_1` must hold n_tok columns of n_embd; with `nw.q8_1_il` 2+ tokens' gate and
+/// up projections run through `native_mmvq_il` (IQ3_S's weights decoded once for all tokens).
 void shared_expert_multi(int n_tok, const float* x, const NativeSharedWeights& nw, const uint16_t* gate_inp_bf16,
                          float* gate, float* up, float* g, float* out, int64_t n_embd, int64_t n_ff, void* stream,
                          void* gate_stream);

src/core/verify.cpp

@@ -479,6 +479,7 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
         nsw.up_type = wsu->native_type; nsw.up_data = wsu->native_data;
         nsw.down_type = wsd->native_type; nsw.down_data = wsd->native_data;
         nsw.q8_1 = xq_;                                   // the main stream's next use is the next layer's mixer
+        nsw.q8_1_il = xil_;
         cudaStream_t ws = cs, gs = cs;
         if (beside) {
             if (cudaStreamWaitEvent(shs_, efork_, 0) != cudaSuccess ||
@@ -528,8 +529,7 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
             else native_quantize_q8_1(x, xq_, (int) n_in, n, cs);
         };
         auto mm = [&](int type, const void* w, float* y, int64_t n_in, int64_t n_out, cudaStream_t ps) {
-            if (n >= 2) native_mmvq_il(type, w, xq_, xil_, y, (int) n_in, (int) n_out, n, ps);
-            else native_mmvq(type, w, xq_, y, (int) n_in, (int) n_out, n, ps);
+            native_mmvq_il(type, w, xq_, xil_, y, (int) n_in, (int) n_out, n, ps);
         };
         auto proj = [&](const WeightRef* w, float* y, int64_t n_in, int64_t n_out, cudaStream_t ps) {
             mm(w->native_type, w->native_data, y, n_in, n_out, ps);
@@ -602,7 +602,8 @@ 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;
                 // beside qkv and the conv: alpha/beta (a few blocks, latency-bound) and the gate projection, joined
-                // before the recurrence
+                // before the recurrence.  qkv's grid is recorded first, so it is dispatched first and the conv runs
+                // beside the gate projection's last blocks.
                 if (!fork(side_, bfork_)) return false;
                 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,
@@ -610,9 +611,9 @@ bool Verifier::record_window(int T, cudaStream_t cs, std::string& err) {
                 if (!mark(side_, bjoin_)) return false;
                 quant(xm, N);
                 if (!fork(shs_, pfork_)) return false;
+                proj(wqkv, qkv + (size_t) tb * C, N, C, cs);
                 proj(wg, z_ + (size_t) tb * ZV, N, ZV, shs_);
                 if (!mark(shs_, pjoin_)) return false;
-                proj(wqkv, qkv + (size_t) tb * C, N, C, cs);
                 gdn_conv_l2_multi(conv, qkv, (const float*) wc->data, hb, (int) C, (int) (2 * HK), EPS, n, cs, tb);
                 if (!wait(bjoin_) || !wait(pjoin_)) return false;
                 // the recurrence from the untouched state over tokens [0, te); outputs only for this group's

src/kernels/cuda/iq_dot.cuh

@@ -245,17 +245,18 @@ struct IQ3SCols {
     }
     template <int NC, class X>
     __device__ __forceinline__ void apply(const X& x, const int& b, float (&out)[NC]) const {
+        int u[8][NC];   // every position's activations first, so that their loads overlap
+#pragma unroll
+        for (int l = 0; l < 8; ++l) x.u(b, l, u[l]);
         int sp[NC], sn[NC];
 #pragma unroll
         for (int c = 0; c < NC; ++c) sp[c] = sn[c] = 0;
 #pragma unroll
         for (int l = 0; l < 8; ++l) {
-            int u[NC];
-            x.u(b, l, u);
 #pragma unroll
             for (int c = 0; c < NC; ++c) {
-                sp[c] = ggml_cuda_dp4a((int) g[l], u[c], sp[c]);
-                sn[c] = ggml_cuda_dp4a((int) gm[l], u[c], sn[c]);
+                sp[c] = ggml_cuda_dp4a((int) g[l], u[l][c], sp[c]);
+                sn[c] = ggml_cuda_dp4a((int) gm[l], u[l][c], sn[c]);
             }
         }
         float d8[NC];

src/kernels/cuda/iq_kernels.cu

@@ -33,43 +33,47 @@ __global__ void __launch_bounds__(128) mmvq_kernel(const uint8_t* __restrict__ w
     }
 }
 
-// mmvq_kernel for NC columns from their interleaved (position-major) copy: a warp takes R rows, decodes each call's
-// weights once and keeps row_dot's lane-strided sum per (row, column), so every value is bitwise mmvq_kernel's.
+// mmvq_kernel for NC columns from their interleaved (position-major) copy, or one from its q8_1 blocks: a warp takes R
+// rows, decodes each call's weights once and keeps row_dot's lane-strided sum per (row, column), so every value is
+// bitwise mmvq_kernel's.
 template <int NC, int R>
 __global__ void __launch_bounds__(128) mmvq_il_iq3_s_kernel(const uint8_t* __restrict__ w, size_t row_bytes,
-                                                            const int* __restrict__ xq, const float* __restrict__ xd,
+                                                            const void* __restrict__ xq, const float* __restrict__ xd,
                                                             float* __restrict__ y, int n_in, int n_out) {
     using F = Fmt<21>;
     const int lane = threadIdx.x & 31;
-    const int row0 = (blockIdx.x * 4 + (threadIdx.x >> 5)) * R;
-    if (row0 >= n_out) return;
-    const Q81Il<NC, true> x{xq, xd, n_in / 32};
+    const auto x = q8_1_cols<NC, true>(xq, xd, n_in / 32);
     const int nb = n_in / F::qk;
-    float acc[R][NC];
+    // a grid-stride loop, though the grid gives each warp one group: native_mmvq_il_kernel's form, which its
+    // rows-a-warp table is measured on
+    for (int grp = blockIdx.x * 4 + (threadIdx.x >> 5); grp * R < n_out; grp += gridDim.x * 4) {
+        const int row0 = grp * R;
+        float acc[R][NC];
 #pragma unroll
-    for (int i = 0; i < R; ++i)
+        for (int i = 0; i < R; ++i)
 #pragma unroll
-        for (int c = 0; c < NC; ++c) acc[i][c] = 0.0f;
-    for (int k = lane; k < nb * F::ipb; k += 32) {
-        const int kbx = k / F::ipb, iqs = F::step * (k % F::ipb);
+            for (int c = 0; c < NC; ++c) acc[i][c] = 0.0f;
+        for (int k = lane; k < nb * F::ipb; k += 32) {
+            const int kbx = k / F::ipb, iqs = F::step * (k % F::ipb);
 #pragma unroll
-        for (int i = 0; i < R; ++i) {
-            const int row = min(row0 + i, n_out - 1);   // a partial last group recomputes its last row
-            IQ3SCols wv;
-            wv.load(w + (size_t) row * row_bytes, kbx, iqs);
-            float o[NC];
-            wv.apply<NC>(x, kbx * (F::qk / 32) + iqs / 2, o);
+            for (int i = 0; i < R; ++i) {
+                const int row = min(row0 + i, n_out - 1);   // a partial last group recomputes its last row
+                IQ3SCols wv;
+                wv.load(w + (size_t) row * row_bytes, kbx, iqs);
+                float o[NC];
+                wv.apply<NC>(x, kbx * (F::qk / 32) + iqs / 2, o);
 #pragma unroll
-            for (int c = 0; c < NC; ++c) acc[i][c] += o[c];
+                for (int c = 0; c < NC; ++c) acc[i][c] += o[c];
+            }
         }
+#pragma unroll
+        for (int i = 0; i < R; ++i)
+#pragma unroll
+            for (int c = 0; c < NC; ++c) {
+                const float s = warp_sum(acc[i][c]);
+                if (lane == 0 && row0 + i < n_out) y[(size_t) c * n_out + row0 + i] = s;
+            }
     }
-#pragma unroll
-    for (int i = 0; i < R; ++i)
-#pragma unroll
-        for (int c = 0; c < NC; ++c) {
-            const float s = warp_sum(acc[i][c]);
-            if (lane == 0 && row0 + i < n_out) y[(size_t) c * n_out + row0 + i] = s;
-        }
 }
 
 // ---------------------------------------------------------------- grouped native experts
@@ -590,21 +594,38 @@ void iq_mmvq(int t, const void* w, const void* x_q8_1, float* y, int n_in, int n
 
 bool iq_mmvq_il_supported(int t) noexcept { return t == 21; }
 
-void iq_mmvq_il(int t, const void* w, const void* x_il, float* y, int n_in, int n_out, int ncols, void* stream) {
-    if (t != 21 || ncols < 2 || ncols > 8 || n_in % 256 != 0) {
+namespace {
+template <int R>
+void launch_iq3_s_il(const uint8_t* w, size_t rb, const void* x_q8_1, const void* x_il, float* y, int n_in,
+                     int n_out, int ncols, cudaStream_t s) {
+    const unsigned blocks = (unsigned) (((n_out + R - 1) / R + 3) / 4);
+    if (ncols == 1) {
+        mmvq_il_iq3_s_kernel<1, R><<<blocks, 128, 0, s>>>(w, rb, x_q8_1, nullptr, y, n_in, n_out);
+        return;
+    }
+    const Q81IlParts parts = q8_1_il_parts(x_il, n_in, ncols);
+    switch (ncols) {
+#define STRATA_IL(N) \
+    case N: mmvq_il_iq3_s_kernel<N, R><<<blocks, 128, 0, s>>>(w, rb, parts.pm, parts.d, y, n_in, n_out); break;
+        STRATA_IL(2) STRATA_IL(3) STRATA_IL(4) STRATA_IL(5) STRATA_IL(6) STRATA_IL(7) STRATA_IL(8)
+#undef STRATA_IL
+    }
+}
+}  // namespace
+
+void iq_mmvq_il(int t, const void* w, const void* x_q8_1, const void* x_il, float* y, int n_in, int n_out, int ncols,
+                void* stream, int rows) {
+    if (t != 21 || ncols < 1 || ncols > 8 || n_in % 256 != 0) {
         std::fprintf(stderr, "iq_mmvq_il: type %d, %d columns are not supported\n", t, ncols);
         std::exit(1);
     }
-    constexpr int R = 2;
-    const Q81IlParts parts = q8_1_il_parts(x_il, n_in, ncols);
     const size_t rb = iq_row_bytes(t, n_in);
     const auto* W = (const uint8_t*) w;
-    const unsigned blocks = (unsigned) (((n_out + R - 1) / R + 3) / 4);
     cudaStream_t s = (cudaStream_t) stream;
-    switch (ncols) {
-#define STRATA_IL(N) case N: mmvq_il_iq3_s_kernel<N, R><<<blocks, 128, 0, s>>>(W, rb, parts.pm, parts.d, y, n_in, n_out); break;
-        STRATA_IL(2) STRATA_IL(3) STRATA_IL(4) STRATA_IL(5) STRATA_IL(6) STRATA_IL(7) STRATA_IL(8)
-#undef STRATA_IL
+    switch (rows) {
+    case 1: launch_iq3_s_il<1>(W, rb, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
+    case 4: launch_iq3_s_il<4>(W, rb, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
+    default: launch_iq3_s_il<2>(W, rb, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
     }
     check("iq_mmvq_il");
 }

src/kernels/cuda/native_mmvq.cu

@@ -30,6 +30,7 @@
 #include <cuda_fp16.h>
 #include <cuda_runtime.h>
 
+#include <algorithm>
 #include <cstdint>
 #include <limits>
 #include <stdexcept>
@@ -1091,8 +1092,8 @@ void launch_multi(const void* weights, const void* x_q8_1, float* y, int n_in, i
 struct IlIQ4XS {
     using F = IQ4XSTraits;
     static constexpr bool PM = true;
-    template <int NC>
-    __device__ static void apply(const F::W& r, const Q81Il<NC, PM>& x, int kby, int iqs, float (&out)[NC]) {
+    template <int NC, class X>
+    __device__ static void apply(const F::W& r, const X& x, int kby, int iqs, float (&out)[NC]) {
         const int b = kby + iqs / 4;
         int sumi[NC];
 #pragma unroll
@@ -1121,8 +1122,8 @@ struct IlIQ4XS {
 struct IlQ4K {
     using F = Q4KTraits;
     static constexpr bool PM = false;
-    template <int NC>
-    __device__ static void apply(const F::W& r, const Q81Il<NC, PM>& x, int kby, int iqs, float (&out)[NC]) {
+    template <int NC, class X>
+    __device__ static void apply(const F::W& r, const X& x, int kby, int iqs, float (&out)[NC]) {
         int u[4][NC];
         float d8[2][NC];
 #pragma unroll
@@ -1155,8 +1156,8 @@ struct IlQ4K {
 struct IlQ5K {
     using F = Q5KTraits;
     static constexpr bool PM = false;
-    template <int NC>
-    __device__ static void apply(const F::W& r, const Q81Il<NC, PM>& x, int kby, int iqs, float (&out)[NC]) {
+    template <int NC, class X>
+    __device__ static void apply(const F::W& r, const X& x, int kby, int iqs, float (&out)[NC]) {
         int u[4][NC];
         float d8[2][NC];
 #pragma unroll
@@ -1193,8 +1194,8 @@ struct IlQ5K {
 struct IlQ6K {
     using F = Q6KTraits;
     static constexpr bool PM = false;
-    template <int NC>
-    __device__ static void apply(const F::W& r, const Q81Il<NC, PM>& x, int kby, int iqs, float (&out)[NC]) {
+    template <int NC, class X>
+    __device__ static void apply(const F::W& r, const X& x, int kby, int iqs, float (&out)[NC]) {
         int u[2][NC];
         float d8[2][NC];
 #pragma unroll
@@ -1231,8 +1232,8 @@ struct IlQ80 {
         r.d0 = w->d;
         return r;
     }
-    template <int NC>
-    __device__ static void apply(const W& r, const Q81Il<NC, PM>& x, int kby, int iqs, float (&out)[NC]) {
+    template <int NC, class X>
+    __device__ static void apply(const W& r, const X& x, int kby, int iqs, float (&out)[NC]) {
         int u0[NC], u1[NC];
         x.u(kby, iqs, u0);
         x.u(kby, iqs + 1, u1);
@@ -1260,67 +1261,85 @@ template <> struct IlLoad<IlQ80> {
 
 template <typename I, int NC, int R>
 __launch_bounds__(4 * WARP)
-__global__ void native_mmvq_il_kernel(const typename I::F::Block* __restrict__ w, const int* __restrict__ xq,
+__global__ void native_mmvq_il_kernel(const typename I::F::Block* __restrict__ w, const void* __restrict__ xq,
                                       const float* __restrict__ xd, float* __restrict__ y, int n_in, int n_out) {
     using F = typename I::F;
     constexpr int SUB = WARP / F::T;   // weight blocks a chunk covers
     const int lane = int(threadIdx.x) & (WARP - 1);
-    const int row0 = (int(blockIdx.x) * 4 + (int(threadIdx.x) >> 5)) * R;
-    if (row0 >= n_out) return;
-    const Q81Il<NC, I::PM> x{xq, xd, n_in / Q8K};
+    const auto x = q8_1_cols<NC, I::PM>(xq, xd, n_in / Q8K);
     const int bpr = n_in / F::DIV;
     const int kqs = F::kqs(lane);
     const int nm = (bpr + SUB - 1) / SUB;
-    float acc[4][R][NC];
+    // A grid-stride loop over the row groups, though the grid gives each warp one: without it ptxas gives several
+    // instances fewer registers, and Q6_K's at 3 columns run 10-50% slower (the rows-a-warp table is measured on it).
+    for (int grp = int(blockIdx.x) * 4 + (int(threadIdx.x) >> 5); grp * R < n_out; grp += int(gridDim.x) * 4) {
+        const int row0 = grp * R;
+        float acc[4][R][NC];
 #pragma unroll
-    for (int v = 0; v < 4; ++v)
+        for (int v = 0; v < 4; ++v)
 #pragma unroll
-        for (int i = 0; i < R; ++i)
+            for (int i = 0; i < R; ++i)
 #pragma unroll
-            for (int c = 0; c < NC; ++c) acc[v][i][c] = 0.0f;
-    for (int m0 = 0; m0 < nm; m0 += 4) {
+                for (int c = 0; c < NC; ++c) acc[v][i][c] = 0.0f;
+        for (int m0 = 0; m0 < nm; m0 += 4) {
 #pragma unroll
-        for (int v = 0; v < 4; ++v) {
-            const int kbx = SUB * (m0 + v) + lane / F::T;
-            if (kbx < bpr) {
+            for (int v = 0; v < 4; ++v) {
+                const int kbx = SUB * (m0 + v) + lane / F::T;
+                if (kbx < bpr) {
 #pragma unroll
-                for (int i = 0; i < R; ++i) {
-                    const int row = min(row0 + i, n_out - 1);   // a partial last group recomputes its last row
-                    const typename IlLoad<I>::W wv = IlLoad<I>::load(w + std::size_t(row) * bpr + kbx, kqs);
-                    float o[NC];
-                    I::template apply<NC>(wv, x, kbx * F::KBY, kqs, o);
+                    for (int i = 0; i < R; ++i) {
+                        const int row = min(row0 + i, n_out - 1);   // a partial last group recomputes its last row
+                        const typename IlLoad<I>::W wv = IlLoad<I>::load(w + std::size_t(row) * bpr + kbx, kqs);
+                        float o[NC];
+                        I::template apply<NC>(wv, x, kbx * F::KBY, kqs, o);
 #pragma unroll
-                    for (int c = 0; c < NC; ++c) acc[v][i][c] += o[c];
+                        for (int c = 0; c < NC; ++c) acc[v][i][c] += o[c];
+                    }
                 }
             }
         }
+#pragma unroll
+        for (int i = 0; i < R; ++i)
+#pragma unroll
+            for (int c = 0; c < NC; ++c) {
+                float s = acc[0][i][c];
+                s += acc[1][i][c];
+                s += acc[2][i][c];
+                s += acc[3][i][c];
+                s = warp_sum(s);
+                if (lane == 0 && row0 + i < n_out) y[std::size_t(c) * n_out + row0 + i] = s;
+            }
     }
-#pragma unroll
-    for (int i = 0; i < R; ++i)
-#pragma unroll
-        for (int c = 0; c < NC; ++c) {
-            float s = acc[0][i][c];
-            s += acc[1][i][c];
-            s += acc[2][i][c];
-            s += acc[3][i][c];
-            s = warp_sum(s);
-            if (lane == 0 && row0 + i < n_out) y[std::size_t(c) * n_out + row0 + i] = s;
-        }
 }
 
 template <typename I, int R>
-void launch_il(const void* weights, const void* x_il, float* y, int n_in, int n_out, int ncols, cudaStream_t s) {
+void launch_il(const void* weights, const void* x_q8_1, const void* x_il, float* y, int n_in, int n_out, int ncols,
+               cudaStream_t s) {
+    const auto* w = static_cast<const typename I::F::Block*>(weights);
+    const unsigned blocks = unsigned(((n_out + R - 1) / R + 3) / 4);
+    if (ncols == 1) {
+        native_mmvq_il_kernel<I, 1, R><<<blocks, 4 * WARP, 0, s>>>(w, x_q8_1, nullptr, y, n_in, n_out);
+        return;
+    }
     const Q81IlParts parts = q8_1_il_parts(x_il, n_in, ncols);
     const int* xq = I::PM ? parts.pm : parts.bm;
     const float* xd = parts.d;
-    const auto* w = static_cast<const typename I::F::Block*>(weights);
-    const unsigned blocks = unsigned(((n_out + R - 1) / R + 3) / 4);
     switch (ncols) {
 #define STRATA_IL_CASE(N) case N: native_mmvq_il_kernel<I, N, R><<<blocks, 4 * WARP, 0, s>>>(w, xq, xd, y, n_in, n_out); break;
         STRATA_IL_CASE(2) STRATA_IL_CASE(3) STRATA_IL_CASE(4) STRATA_IL_CASE(5) STRATA_IL_CASE(6) STRATA_IL_CASE(7)
         STRATA_IL_CASE(8)
 #undef STRATA_IL_CASE
-        default: throw std::invalid_argument("native_mmvq_il requires 2 <= ncols <= 8");
+        default: throw std::invalid_argument("native_mmvq_il requires 1 <= ncols <= 8");
+    }
+}
+
+template <typename I>
+void launch_il_rows(int r, const void* weights, const void* x_q8_1, const void* x_il, float* y, int n_in, int n_out,
+                    int ncols, cudaStream_t s) {
+    switch (r) {
+    case 1: launch_il<I, 1>(weights, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
+    case 4: launch_il<I, 4>(weights, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
+    default: launch_il<I, 2>(weights, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
     }
 }
 
@@ -1781,37 +1800,60 @@ void native_quantize_q8_1_il(const float* x, void* x_q8_1, void* x_il, int n_in,
     launch_check();
 }
 
+namespace {
+// Rows a warp at 1-4 columns ([ncols - 1]) by the row count, in classes around the model's shapes: below 2048 (the
+// shared expert's 640, the keys' 512), 4096 (the output projections' 2560), 8192 (the gate projection's 6144), 12288
+// (qkv's 10240), and more (the queries' 12288, the head's); 0: native_mmvq's kernels, which spread a row over four
+// warps (IQ3_S's decode the weights once a column).  The fastest in the verify window on the RTX 3090
+// (scratch/mm1_bench.cu alone, CUPTI traces beside the other branches); 5-8 columns take 2 rows from 2048.
+struct IlRows { int type; uint8_t r[4][5]; };
+constexpr IlRows kIlRows[] = {
+    {23, {{0, 0, 0, 0, 0}, {0, 2, 2, 2, 2}, {0, 2, 2, 2, 2}, {0, 2, 2, 1, 2}}},   // IQ4_XS
+    {12, {{0, 0, 0, 0, 0}, {0, 2, 2, 2, 2}, {0, 2, 2, 2, 2}, {0, 2, 2, 2, 2}}},   // Q4_K
+    {13, {{0, 0, 0, 0, 0}, {0, 2, 1, 1, 2}, {0, 2, 2, 2, 2}, {0, 2, 2, 2, 2}}},   // Q5_K
+    {14, {{0, 0, 0, 0, 0}, {0, 1, 1, 1, 1}, {0, 1, 1, 1, 1}, {0, 2, 2, 1, 2}}},   // Q6_K
+    {8, {{0, 0, 0, 0, 0}, {0, 0, 4, 2, 1}, {0, 2, 2, 2, 2}, {0, 2, 2, 2, 2}}},    // Q8_0
+    {21, {{1, 1, 1, 2, 2}, {1, 2, 4, 2, 2}, {1, 1, 1, 1, 1}, {1, 2, 2, 1, 1}}},   // IQ3_S
+};
+int il_rows(int type, int ncols, int n_out) {
+    const int cls = n_out < 2048 ? 0 : n_out < 4096 ? 1 : n_out < 8192 ? 2 : n_out < 12288 ? 3 : 4;
+    for (const IlRows& e : kIlRows)
+        if (e.type == type) return ncols <= 4 ? e.r[ncols - 1][cls] : cls ? 2 : 0;
+    return 0;
+}
+
+int g_tune_rows = 0;         // native_mmvq_il_tune (tests, benchmarks): rows a warp for every shape (0: the table),
+bool g_tune_one = false;     // and 1 column through these kernels too
+}  // namespace
+void native_mmvq_il_tune(int rows, bool one) {
+    g_tune_rows = rows;
+    g_tune_one = one;
+}
+
 void native_mmvq_il(int ggml_type, const void* weights, const void* x_q8_1, const void* x_il, float* y, int n_in,
                     int n_out, int ncols, void* stream) {
-    // below ~2048 rows the warps are too few for one row each: the multi-column kernel spreads a row over four (and
-    // for Q8_0 2 columns of fewer than 4096 rows)
     const bool il = ggml_type == 23 || ggml_type == 12 || ggml_type == 13 || ggml_type == 14 || ggml_type == 8;
-    if (ncols < 2 || n_out < 2048 || (ggml_type == 8 && ncols == 2 && n_out < 4096) ||
-        (!il && !iq_mmvq_il_supported(ggml_type))) {
+    const int r = g_tune_rows && (ncols > 1 || g_tune_one) ? g_tune_rows : il_rows(ggml_type, ncols, n_out);
+    if (r == 0 || (!il && !iq_mmvq_il_supported(ggml_type))) {
         native_mmvq(ggml_type, weights, x_q8_1, y, n_in, n_out, ncols, stream);
         return;
     }
     if (!il) {
-        iq_mmvq_il(ggml_type, weights, x_il, y, n_in, n_out, ncols, stream);
+        iq_mmvq_il(ggml_type, weights, x_q8_1, x_il, y, n_in, n_out, ncols, stream, r);
         return;
     }
     validate_shape(n_in, ncols, ggml_type == 8 ? Q8K : QK);
     validate_pointer(weights);
-    validate_pointer(x_il);
+    validate_pointer(ncols == 1 ? x_q8_1 : x_il);
     validate_pointer(y);
     validate_stream(stream);
     const auto s = static_cast<cudaStream_t>(stream);
-    // rows a warp: 2, and 4 for Q4_K and Q6_K from 8192 rows (the faster on the model files' shapes, mmvq_il_parity
-    // --bench)
-    const bool wide = n_out >= 8192;
     switch (ggml_type) {
-    case 23: launch_il<IlIQ4XS, 2>(weights, x_il, y, n_in, n_out, ncols, s); break;
-    case 12: wide ? launch_il<IlQ4K, 4>(weights, x_il, y, n_in, n_out, ncols, s)
-                  : launch_il<IlQ4K, 2>(weights, x_il, y, n_in, n_out, ncols, s); break;
-    case 13: launch_il<IlQ5K, 2>(weights, x_il, y, n_in, n_out, ncols, s); break;
-    case 14: wide ? launch_il<IlQ6K, 4>(weights, x_il, y, n_in, n_out, ncols, s)
-                  : launch_il<IlQ6K, 2>(weights, x_il, y, n_in, n_out, ncols, s); break;
-    case 8: launch_il<IlQ80, 2>(weights, x_il, y, n_in, n_out, ncols, s); break;
+    case 23: launch_il_rows<IlIQ4XS>(r, weights, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
+    case 12: launch_il_rows<IlQ4K>(r, weights, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
+    case 13: launch_il_rows<IlQ5K>(r, weights, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
+    case 14: launch_il_rows<IlQ6K>(r, weights, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
+    case 8: launch_il_rows<IlQ80>(r, weights, x_q8_1, x_il, y, n_in, n_out, ncols, s); break;
     }
     launch_check();
 }

src/kernels/cuda/q8_1_il.cuh

@@ -4,9 +4,11 @@
 
 #include "strata/kernels/native_mmvq.hpp"
 
+#include <cuda_fp16.h>
 #include <cuda_runtime.h>
 
 #include <cstddef>
+#include <cstdint>
 
 namespace strata::kernels {
 namespace {
@@ -62,6 +64,25 @@ template <int NC, bool PM> struct Q81Il {
     }
 };
 
+// One column's plain q8_1 blocks (36 bytes: the half2 scale and sum, then 32 codes) behind Q81Il's interface
+struct Q81One {
+    const uint8_t* __restrict__ x;
+    __device__ __forceinline__ void u(int b, int p, int (&o)[1]) const {
+        o[0] = reinterpret_cast<const int*>(x + std::size_t(b) * 36 + 4)[p];
+    }
+    __device__ __forceinline__ void scales(int b, float (&o)[1]) const {
+        o[0] = __low2float(*reinterpret_cast<const half2*>(x + std::size_t(b) * 36));
+    }
+};
+
+// The activations a kernel for NC columns reads: one column's q8_1 blocks, 2+ columns' interleaved copy (nb blocks a
+// column)
+template <int NC, bool PM>
+__device__ __forceinline__ auto q8_1_cols(const void* xq, const float* xd, int nb) {
+    if constexpr (NC == 1) return Q81One{static_cast<const uint8_t*>(xq)};
+    else return Q81Il<NC, PM>{static_cast<const int*>(xq), xd, nb};
+}
+
 // The copy's parts inside the buffer native_q8_1_il_bytes sizes: bm, pm (int32) and d (float).
 struct Q81IlParts {
     const int* bm;

src/kernels/cuda/shared_expert.cu

@@ -176,9 +176,16 @@ void shared_expert_multi(int n_tok, const float* x, const NativeSharedWeights& n
                                     "streams");
     cudaStream_t cs = (cudaStream_t) stream;
     bf16_gemv_fp32_mmvf_multi(x, gate_inp_bf16, g, n_embd, 1, n_tok, gate_stream);
-    native_quantize_q8_1(x, nw.q8_1, (int) n_embd, n_tok, stream);
-    native_mmvq(nw.gate_type, nw.gate_data, nw.q8_1, gate, (int) n_embd, (int) n_ff, n_tok, stream);
-    native_mmvq(nw.up_type, nw.up_data, nw.q8_1, up, (int) n_embd, (int) n_ff, n_tok, stream);
+    // native_mmvq_il's kernels where they are the faster (1 token, or 2+ with the interleaved copy)
+    const bool il = n_tok == 1 || nw.q8_1_il;
+    if (n_tok >= 2 && nw.q8_1_il) native_quantize_q8_1_il(x, nw.q8_1, nw.q8_1_il, (int) n_embd, n_tok, stream);
+    else native_quantize_q8_1(x, nw.q8_1, (int) n_embd, n_tok, stream);
+    auto proj = [&](int type, const void* w, float* y) {
+        if (il) native_mmvq_il(type, w, nw.q8_1, nw.q8_1_il, y, (int) n_embd, (int) n_ff, n_tok, stream);
+        else native_mmvq(type, w, nw.q8_1, y, (int) n_embd, (int) n_ff, n_tok, stream);
+    };
+    proj(nw.gate_type, nw.gate_data, gate);
+    proj(nw.up_type, nw.up_data, up);
     const int n = (int) (n_ff * n_tok);
     native_swiglu_kernel<<<(unsigned) ((n + THREADS - 1) / THREADS), THREADS, 0, cs>>>(gate, up, gate, n);
     native_quantize_q8_1(gate, nw.q8_1, (int) n_ff, n_tok, stream);

src/kernels/mmvq_il_parity.cu

@@ -1,6 +1,7 @@
-// src/kernels/mmvq_il_parity.cu - native_mmvq_il (2..8 columns from the interleaved q8_1 copy) against native_mmvq,
-// bitwise, on synthetic matrices of the model files' dense formats and shapes; the interleaving quantizer's q8_1 blocks
-// against native_quantize_q8_1's.  --bench times both inside CUDA graphs (weights cycled past the L2).
+// src/kernels/mmvq_il_parity.cu - native_mmvq_il (1 column from its q8_1 blocks, 2..8 from the interleaved copy)
+// against native_mmvq, bitwise, on synthetic matrices of the model files' dense formats and shapes, with every rows a
+// warp; the interleaving quantizer's q8_1 blocks against native_quantize_q8_1's.  --bench times both inside CUDA graphs
+// (weights cycled past the L2).
 #include "strata/kernels/native_mmvq.hpp"
 
 #include <cuda_runtime.h>
@@ -78,6 +79,7 @@ std::vector<T> fetch(const void* d, size_t n) {
     return h;
 }
 
+// every column count 1..8 (1: the one-column kernels), with the table's rows a warp and each of 1, 2, 4
 int check_shape(const Fmt& f, const Act& a, int n_out) {
     const int n_in = a.n_in;
     void* w = dev_weights(f, n_in, n_out);
@@ -85,25 +87,33 @@ int check_shape(const Fmt& f, const Act& a, int n_out) {
     CK(cudaMalloc(&y0, (size_t) 8 * n_out * 4));
     CK(cudaMalloc(&y1, (size_t) 8 * n_out * 4));
     int bad = 0;
-    for (int nc = 2; nc <= 8; ++nc) {
+    for (int nc = 1; nc <= 8; ++nc) {
         native_quantize_q8_1(a.x, a.q_ref, n_in, nc, g_s);
-        native_quantize_q8_1_il(a.x, a.q, a.il, n_in, nc, g_s);
+        if (nc >= 2) native_quantize_q8_1_il(a.x, a.q, a.il, n_in, nc, g_s);
+        else native_quantize_q8_1(a.x, a.q, n_in, nc, g_s);
         CK(cudaMemsetAsync(y0, 0xff, (size_t) 8 * n_out * 4, g_s));
-        CK(cudaMemsetAsync(y1, 0xee, (size_t) 8 * n_out * 4, g_s));
         native_mmvq(f.type, w, a.q_ref, y0, n_in, n_out, nc, g_s);
-        native_mmvq_il(f.type, w, a.q, a.il, y1, n_in, n_out, nc, g_s);
         CK(cudaStreamSynchronize(g_s));
         const size_t qb = native_q8_1_bytes(n_in, nc);
         if (fetch<uint8_t>(a.q_ref, qb) != fetch<uint8_t>(a.q, qb)) {
             ++bad;
             std::printf("  %s %d->%d, %d columns: the q8_1 blocks differ\n", f.name, n_in, n_out, nc);
         }
-        const auto r0 = fetch<uint32_t>(y0, (size_t) nc * n_out), r1 = fetch<uint32_t>(y1, (size_t) nc * n_out);
-        size_t diff = 0;
-        for (size_t i = 0; i < r0.size(); ++i) diff += r0[i] != r1[i];
-        if (diff) {
-            ++bad;
-            std::printf("  %s %d->%d, %d columns: %zu of %zu values differ\n", f.name, n_in, n_out, nc, diff, r0.size());
+        const auto r0 = fetch<uint32_t>(y0, (size_t) nc * n_out);
+        for (int rows : {0, 1, 2, 4}) {
+            CK(cudaMemsetAsync(y1, 0xee, (size_t) 8 * n_out * 4, g_s));
+            native_mmvq_il_tune(rows, nc == 1);
+            native_mmvq_il(f.type, w, a.q, a.il, y1, n_in, n_out, nc, g_s);
+            native_mmvq_il_tune(0, false);
+            CK(cudaStreamSynchronize(g_s));
+            const auto r1 = fetch<uint32_t>(y1, (size_t) nc * n_out);
+            size_t diff = 0;
+            for (size_t i = 0; i < r0.size(); ++i) diff += r0[i] != r1[i];
+            if (diff) {
+                ++bad;
+                std::printf("  %s %d->%d, %d columns, rows %d: %zu of %zu values differ\n", f.name, n_in, n_out, nc,
+                            rows, diff, r0.size());
+            }
         }
     }
     cudaFree(w);
@@ -112,6 +122,13 @@ int check_shape(const Fmt& f, const Act& a, int n_out) {
     return bad;
 }
 
+__global__ void spin_kernel(long long cycles) {
+    const long long t0 = clock64();
+    while (clock64() - t0 < cycles) {}
+}
+
+// The median of 20 launches, each after ~3 ms of a one-thread spin, as the engine's load: the GPU active at little
+// power (a continuous load reaches the 3090's 260 W cap and lowers the clocks, an idle GPU lowers them too).
 double time_calls(const std::function<void(int)>& fn, int reps = 48) {
     cudaGraph_t g;
     CK(cudaStreamBeginCapture(g_s, cudaStreamCaptureModeThreadLocal));
@@ -120,25 +137,25 @@ double time_calls(const std::function<void(int)>& fn, int reps = 48) {
     cudaGraphExec_t ex;
     CK(cudaGraphInstantiate(&ex, g, 0));
     cudaGraphDestroy(g);
-    const auto t0 = std::chrono::steady_clock::now();
-    while (std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count() < 0.3) {   // clocks up
-        CK(cudaGraphLaunch(ex, g_s));
-        CK(cudaStreamSynchronize(g_s));
-    }
     cudaEvent_t a, b;
     CK(cudaEventCreate(&a));
     CK(cudaEventCreate(&b));
-    const int iters = 20;
-    CK(cudaEventRecord(a, g_s));
-    for (int i = 0; i < iters; ++i) CK(cudaGraphLaunch(ex, g_s));
-    CK(cudaEventRecord(b, g_s));
-    CK(cudaEventSynchronize(b));
-    float ms = 0;
-    CK(cudaEventElapsedTime(&ms, a, b));
+    std::vector<double> ts;
+    for (int k = 0; k < 24; ++k) {
+        spin_kernel<<<1, 1, 0, g_s>>>(5'500'000);
+        CK(cudaEventRecord(a, g_s));
+        CK(cudaGraphLaunch(ex, g_s));
+        CK(cudaEventRecord(b, g_s));
+        CK(cudaEventSynchronize(b));
+        float ms = 0;
+        CK(cudaEventElapsedTime(&ms, a, b));
+        if (k >= 4) ts.push_back(ms * 1000.0 / reps);
+    }
+    std::sort(ts.begin(), ts.end());
     cudaGraphExecDestroy(ex);
     cudaEventDestroy(a);
     cudaEventDestroy(b);
-    return ms * 1000.0 / iters / reps;
+    return ts[ts.size() / 2];
 }
 
 void bench_shape(const Fmt& f, const Act& a, int n_out) {
@@ -150,7 +167,7 @@ void bench_shape(const Fmt& f, const Act& a, int n_out) {
     float* y;
     CK(cudaMalloc(&y, (size_t) 8 * n_out * 4));
     std::printf("%-6s %5d -> %5d (%5.2f MB), us (GB/s):", f.name, n_in, n_out, mb);
-    for (int nc : {2, 3, 4, 8}) {
+    for (int nc : {2, 3, 4}) {
         native_quantize_q8_1(a.x, a.q_ref, n_in, nc, g_s);
         native_quantize_q8_1_il(a.x, a.q, a.il, n_in, nc, g_s);
         const double t0 = time_calls([&](int r) { native_mmvq(f.type, w[(size_t) (r % copies)], a.q_ref, y, n_in, n_out, nc, g_s); });
@@ -172,16 +189,19 @@ int main(int argc, char** argv) {
     Act a2560 = make_act(2560), a6144 = make_act(6144), a640 = make_act(640);
     if (bench) {
         for (const Fmt& f : kFmts) {
+            bench_shape(f, a2560, 12288);
             bench_shape(f, a2560, 10240);
             bench_shape(f, a2560, 6144);
             bench_shape(f, a6144, 2560);
+            bench_shape(f, a2560, 640);
         }
         return 0;
     }
     int bad = 0, cases = 0;
     for (const Fmt& f : kFmts) {
         const std::pair<Act*, int> shapes[] = {{&a2560, 10240}, {&a2560, 6144}, {&a6144, 2560}, {&a2560, 12288},
-                                               {&a2560, 2050}, {&a2560, 4097}, {&a2560, 8192}, {&a2560, 512}};
+                                               {&a2560, 2050}, {&a2560, 4097}, {&a2560, 8192}, {&a2560, 512},
+                                               {&a2560, 640}};
         for (const auto& sh : shapes) {
             bad += check_shape(f, *sh.first, sh.second);
             ++cases;
@@ -191,6 +211,7 @@ int main(int argc, char** argv) {
             ++cases;
         }
     }
-    std::printf("mmvq_il_parity: %d shape/format cases x 2..8 columns: %s\n", cases, bad ? "FAILED" : "bitwise equal");
+    std::printf("mmvq_il_parity: %d shape/format cases x 1..8 columns x rows a warp: %s\n", cases,
+                bad ? "FAILED" : "bitwise equal");
     return bad ? 1 : 0;
 }