sycl: QSA block scores as oneMKL GEMM tiles from 8K cells (selection 3x faster at 40K, prompt +21% at 40K, +35% at 128K)

The prompt path scores every pooled key block against every query's indexer heads; on SYCL that is the warp
kernel (the CUDA build's tensor-core path is not ported), 8% of a 40K prompt on a B70 and growing with the
context. As a dense product it is one oneMKL FP32 GEMM per tile (pooled keys x the batch's 4 x 256 indexer
queries) and a small kernel that applies the relu per head, sums in head order, and gives block n_bid the dead
key and its +1e9 rule (qsa_block_scores_reduce, the warp kernel's semantics). The tile is carved from the
attention scratch, free until that layer's attention on the same queue, so no VRAM is added. Below ~8K cells
oneMKL's launch cost loses (-2% at 2K, even at 8K), so by default only a batch reaching 2,048 blocks takes it;
STRATA_SELECT_GEMM=0 never, =1 always.

FP32-level, not bitwise (as the CUDA 3xTF32 path). Arc Pro B70, `--pcie-frac 0 --adapt-every 0
STRATA_PREFILL_MMQ=1`, with the per-cell attention, old kernels -> new: Coder IQ1_M 40K 1,562 -> 1,885 tok/s,
128K 1,336 -> 1,798 (the 256 greedy tokens identical); IQ2_XS 40K 1,373 -> 1,613. First-token KL against the
old kernels at 40K 1.3e-11 nats, top-1 equal. docs/INTEL.md has the table.

docs/INTEL.md

@@ -979,3 +979,27 @@ medians, output tokens identical in every arm):
 
 The B70's int8 expert kernels also build and pass the parity test as a SPIR-V JIT build and as an AOT build for `bmg-g21` (compile only: no
 G21 card here); an AOT build for another device (Alchemist, Xe-HPC) gets stubs (`STRATA_I8_NO_ESIMD`).
+
+## Arc Pro B70: QSA prompt attention per cell and block scores as GEMM tiles (2026-10-19)
+
+Two prompt-path kernels that grow with the context. On a 40K prompt (Coder IQ1_M, 256K config) the QSA attention was 27.6%
+of the GPU timeline (7.3 s) and the block selection 8.0% (2.1 s).
+
+- **Attention scores one work-item per cell** (`qsa_prompt_attn_percell`, prompt path only; `STRATA_ATTN_PERCELL=0` off): a
+  work-item holds one cell's key row and all 12 heads' dot products instead of a sub-group splitting the 256 dimensions with a
+  shuffle tree. 2.0-2.2x the warp kernel in `qsa_prompt_attn_parity` (int8 / FP16 pools; 1.2-1.4x on Q4_0).
+- **Block scores as oneMKL FP32 GEMM tiles** (pooled keys x the batch's indexer queries, then a relu-sum kernel;
+  `qsa_block_scores_reduce`). The tile lives in the attention scratch, so it costs no VRAM. It loses to the warp kernel below
+  ~8K cells (oneMKL's launch: -2% at 2K, even at 8K), so it runs only for a batch reaching 2,048 blocks; `STRATA_SELECT_GEMM=0`
+  never, `=1` always. Selection at 40K 2.1 -> 0.7 s.
+
+Both are FP32 in another summation order, not bitwise. First-token logits against the old kernels (whole vocabulary): KL
+3.9e-10 nats at 8K and 1.3e-11 at 40K, top-1 equal (upstream's int8 against FP16 prompt experts, the same prompts: 8.4e-10 and
+5.1e-12); the 256 greedy tokens after a 128K prompt are identical.
+
+`--pcie-frac 0 --adapt-every 0 STRATA_PREFILL_MMQ=1`, 256 greedy tokens, interleaved runs, prefill tok/s (old -> new):
+
+| model | 2,185 | 40,000 | 128,000 |
+|---|---|---|---|
+| Coder IQ1_M, 256K config | 1,329 -> 1,468 (+10%) | 1,562 -> 1,885 (+21%) | 1,336 -> 1,798 (+35%) |
+| Flash-Next IQ2_XS, 128K config | | 1,373 -> 1,613 (+18%) | |

sycl/src/kernels/cuda/qsa_select.dp.cpp

@@ -1566,6 +1566,36 @@ void block_topk_cluster_kernel(const float *__restrict__ scores,
 #endif  // !__HIPCC__
 }  // namespace
 
+// SYCL port: the prompt path's block scores as oneMKL GEMM tiles (prefill.cpp): S[q * IDX_HEADS + h][b - b0] is
+// pooled[b] . q_idx[q][h] for the blocks [b0, b0 + nb); this turns them into scores[q][b] - relu per indexer head,
+// summed in head order - and gives block n_bid the dead key with block_scores_kernel's +1e9 rule. Blocks past a
+// query's n_bid are not written, as in qsa_block_scores. FP32, another summation order than the warp kernel: not
+// bitwise (relative error 2-4e-7; the same selection but for near-ties).
+void qsa_block_scores_reduce(const float* S, int64_t ld, int64_t b0, int64_t nb, const float* dead, const float* q_idx,
+                             const int32_t* steps, int64_t nq, int64_t max_blocks, float* scores, void* stream) {
+    if (nq <= 0 || nb <= 0) return;
+    strata::q_of(stream)->parallel_for(sycl::range<2>((size_t) nq, (size_t) nb), [=](sycl::id<2> it) {
+        const int64_t qi = (int64_t) it[0], b = b0 + (int64_t) it[1];
+        const int64_t n_kv = steps[qi * kStepCount + kStepNKv], n_bid = steps[qi * kStepCount + kStepNBid];
+        if (b > n_bid || b >= max_blocks) return;
+        float score = 0.0f;
+        if (b == n_bid) {
+            for (int h = 0; h < IDX_HEADS; ++h) {
+                float d = 0.0f;
+                for (int i = 0; i < IDX_DIM; ++i) d += dead[i] * q_idx[(qi * IDX_HEADS + h) * IDX_DIM + i];
+                score += d > 0.0f ? d : 0.0f;
+            }
+            if (n_kv % R != 0) score += 1e9f;
+        } else {
+            for (int h = 0; h < IDX_HEADS; ++h) {
+                const float d = S[(size_t) (qi * IDX_HEADS + h) * (size_t) ld + (size_t) (b - b0)];
+                score += d > 0.0f ? d : 0.0f;
+            }
+        }
+        scores[(size_t) qi * (size_t) max_blocks + (size_t) b] = score;
+    });
+}
+
 void qsa_block_scores(const float* pooled, const float* dead, const float* q_idx, const int32_t* steps, int64_t nq,
                       int64_t max_blocks, const QsaShapes& s, float* scores, void* stream, int64_t active_blocks) {
     if (nq <= 0) return;

sycl/src/prefill/prefill.cpp

@@ -40,7 +40,10 @@
 #include "strata/prefill/moe_mmq.hpp"
 #include "strata/core/peer_experts.hpp"
 #include "strata/prefill/kernels.hpp"
-namespace strata::kernels {   // SYCL port, kept out of the shared headers: qsa_decode_attn.dp.cpp
+#include <oneapi/mkl.hpp>   // SYCL port: the QSA block scores as GEMM tiles
+namespace strata::kernels {   // SYCL port, kept out of the shared headers: qsa_select.dp.cpp, qsa_decode_attn.dp.cpp
+void qsa_block_scores_reduce(const float* S, int64_t ld, int64_t b0, int64_t nb, const float* dead, const float* q_idx,
+                             const int32_t* steps, int64_t nq, int64_t max_blocks, float* scores, void* stream);
 void qsa_prompt_attn_percell(const float* q, const QsaAttnPools& pools, const int32_t* ids, const int32_t* steps,
                              int64_t cap, const QsaShapes& s, float* scratch, float* attn, int64_t n_q, void* stream);
 }
@@ -3274,9 +3277,41 @@ bool Prefill::run_impl(const int64_t *tokens, int64_t n, int64_t pos0,
                                                                                 strata::kernels::kStepNBid)] + 1;
                         // the scores on tensor cores (3xTF32: FP32-level, not bitwise); STRATA_SELECT_OLD=1: the warp kernel
                         static const bool old_sel = std::getenv("STRATA_SELECT_OLD") != nullptr;
-                        if (old_sel || !strata::kernels::qsa_block_scores_tc(st.idx_pooled, st.idx_dead, m.q_idx + t0 * 512,
+                        // SYCL port: the scores as oneMKL FP32 GEMM tiles (pooled keys x the batch's indexer queries)
+                        // and a relu-sum (qsa_block_scores_reduce), FP32-level like the CUDA build's 3xTF32 path. The
+                        // tile lives in the attention scratch, which is free until this layer's attention (same
+                        // queue). It pays from ~8K cells on a B70 (launch cost below: -2% at 2K, even at 8K, select
+                        // 3x faster at 40K), so by default only for a batch reaching 2,048 blocks.
+                        // STRATA_SELECT_GEMM=0: the warp kernel always, =1: the GEMM always.
+                        static const int sel_gemm_mode = [] { const char* v = std::getenv("STRATA_SELECT_GEMM"); return v ? std::atoi(v) : -1; }();
+                        const bool sel_gemm = sel_gemm_mode > 0 || (sel_gemm_mode < 0 && active >= 2048);
+                        const int64_t sel_tile = std::min<int64_t>(8192, (int64_t) ((uint64_t) m.attn_batch *
+                            strata::kernels::qsa_decode_attn_scratch_floats(m.cap, s) / (uint64_t) (nb * 4)));
+                        bool sel_done = false;
+                        if (sel_gemm && !old_sel && m.attn_scratch != nullptr && s.idx_dim == 128 && s.idx_n_head == 4 &&
+                            sel_tile >= 256) {
+                            try {
+                                for (int64_t b0 = 0; b0 < active; b0 += sel_tile) {
+                                    const int64_t nbk = std::min<int64_t>(sel_tile, active - b0);
+                                    oneapi::mkl::blas::column_major::gemm(*m.cs, oneapi::mkl::transpose::trans,
+                                                                          oneapi::mkl::transpose::nontrans, nbk, nb * 4, 128,
+                                                                          1.0f, st.idx_pooled + b0 * 128, 128,
+                                                                          m.q_idx + t0 * 512, 128, 0.0f, m.attn_scratch,
+                                                                          sel_tile);
+                                    strata::kernels::qsa_block_scores_reduce(m.attn_scratch, sel_tile, b0, nbk, st.idx_dead,
+                                                                             m.q_idx + t0 * 512, steps0, nb, m.max_blocks,
+                                                                             m.sel_scores, m.cs);
+                                }
+                                sel_done = true;
+                            } catch (const std::exception& ex) {
+                                static bool said = false;
+                                if (!said) std::fprintf(stderr, "prefill: the GEMM block scores failed (%s); the warp kernel instead\n", ex.what());
+                                said = true;
+                            }
+                        }
+                        if (!sel_done && (old_sel || !strata::kernels::qsa_block_scores_tc(st.idx_pooled, st.idx_dead, m.q_idx + t0 * 512,
                                                                              steps0, nb, m.max_blocks, s, m.sel_scores,
-                                                                             m.cs, active))
+                                                                             m.cs, active)))
                             strata::kernels::qsa_block_scores(st.idx_pooled, st.idx_dead, m.q_idx + t0 * 512, steps0, nb,
                                                               m.max_blocks, s, m.sel_scores, m.cs, active);
                         strata::kernels::qsa_block_topk(m.sel_scores, steps0, nb, m.max_blocks, m.cap, s,