Virtual CSR over the DR cube; aggregation as BF16 GEMM (ndarray tile + ndarray CBLAS); benchmark of every arm

.claude/board/entries/2026-11-08-dr-virtual-csr-bf16-gemm.md

@@ -0,0 +1,112 @@
+# 2026-11-08 — Virtual CSR over the DR cube; aggregation as BF16 GEMM (tile, ndarray CBLAS)
+
+**Status:** MEASURED (lab). This **fails the ms gate**: every aggregation row
+is in milliseconds. Crate: `crates/dr-isa-stream-lab` (`vcsr.rs`, `spmm.rs`,
+`bin/bench_vcsr.rs`). Stacked on the DR × ISA stream lab entry of the same
+day. Depends on AdaWorldAPI/ndarray#414: the strided GEMM entry.
+
+## What it is
+
+**A virtual CSR.** Adjacency lives in the lanes' own 256 × 256 cube
+coordinates.
+
+- Each row tile lists its nonempty column tiles.
+- Each (row tile, column tile) pair is one 256 × 256 bitmask block.
+- There is no `row_ptr` or `col_idx` array. Row `i`'s neighbours are the set
+  bits of block row `slot(i)` across the blocks of `tile(i)`.
+- `ExplicitCsr` is the materialised form of the same edges. It is only an
+  oracle, and the row-by-row test shows the two agree.
+
+**Neighbour aggregation (plus-times).** `Y[p][i] = Σ_{j ∈ row(i)} lane_p[j]`
+for 16 bound lanes, computed per block as BF16 GEMMs:
+
+- **ndarray tile GEMM** (AMX TDPBF16PS on this host). A is read in place from
+  the slab with stride = pitch. B is staged per 16-row group in VNNI layout,
+  and empty groups are skipped.
+- **ndarray's own** `backend::cblas_gemm_bf16bf16f32`: the CBLAS-compatible
+  entry (AMX, then AVX-512 BF16, then scalar), with no external library. It
+  takes dense `&[u16]` and no `lda`, so A is staged; one call per block. (An
+  earlier revision of this lab used upstream `mkl_rt` instead. That arm was
+  removed; the comparator is ndarray's own GEMM.)
+- **Staged A**, for comparison.
+
+All three are exact against the CSR oracle at n = 256, 513 (a partial last
+tile) and 4096 for the tile paths. The tile paths are also exact at 65,536,
+and ndarray's CBLAS entry is exact at 513 and 65,536.
+
+**Not claimed:** blasgraph's `CsrStorage` holds `BitVec` values under HDR
+semirings (XorBundle, HammingMin, …). A floating-point GEMM does not realise
+those, and they were not wired: depending on `lance-graph` pulls in DataFusion
+and Lance unconditionally.
+
+## Measured
+
+See `crates/dr-isa-stream-lab/BENCH.md`, which has every arm of both labs
+in one table, from one run. For the graph rows, ndarray's CBLAS entry costs
+581 ms, 4.7 s and 51 ms on the three graphs, because it repacks both
+operands on every call (see BENCH.md, Reading).
+
+## What the numbers say
+
+1. **The cost of a block GEMM is set by block area, not by edges.** Each block
+   costs 16 × 256 × 256 MACs regardless of content.
+   - At ~101 edges per block (0.15% fill), the GEMM loses 2× to a scalar CSR
+     loop over the same lanes.
+   - At 25% fill (community structure), it wins 6.4×.
+   - This is the same crossover the WHERE-mask GEMM showed from the other
+     side: there B was ~50% dense, and the GEMM was the right tool.
+2. **The block layout is only as sparse as the local order is local.** The
+   same degree with one long-range edge per identity grows from 5,172 to
+   41,916 blocks (8.1×) and 330 MB of bitmask. A virtual CSR therefore wants a
+   locality-preserving G2L order for graph workloads. Nothing in DR provides
+   that order yet; it is a property of how local ordinals are assigned.
+3. **Reading A in place does not matter here.** Direct and staged A are within
+   5–10%: per call, B dominates, not A.
+4. **Per call, the community case costs ~3 µs** against the ~0.45 µs measured
+   for a bare K = 256 GEMM. The likely cause is restaging B for every call
+   (1,024 VNNI writes plus clears), though the adjacency is constant per graph
+   and could be packed once. **Not measured separately.**
+5. **ndarray's CBLAS entry is 3–45× slower than the strided tile path.** It
+   repacks A and B on every call: a prepacked or strided entry is the gap.
+
+## Production crate: not now (DECISION)
+
+**DECISION:** no production crate for the virtual CSR or the graph
+aggregation.
+
+**BASIS (measured, BENCH.md):**
+
+- Every graph-aggregation arm misses the gate, and the best one loses to a
+  scalar CSR loop on the sparse graph.
+- Even the dense community case is bounded below by its own work:
+  16 × 65,536 × 256 MACs ≈ 268M. At the measured AMX rate (~43 ns per
+  8,192-MAC block) that is ≈ 1.4 ms before any staging. So the batch cannot
+  reach µs at this size, whichever arm runs it.
+- The one arm that meets the gate is `SUM(lane) WHERE mask` over the whole
+  universe: 172 µs per 16 lanes × 16 masks with A in place. That is an
+  aggregation terminal over masks. If it is promoted, its home is the
+  Quack / mask-risc terminal set, behind DT64 binding, not a new crate.
+
+**SCOPE:** `crates/dr-isa-stream-lab` stays a lab.
+
+**REVISIT WHEN** one of these is measured:
+
+- each block's B packed once per graph version;
+- a per-block dense/sparse choice;
+- a locality-preserving G2L order;
+- a prepacked or strided entry in ndarray's CBLAS GEMM;
+
+and a graph arm then beats the scalar CSR loop on a sparse graph, or the
+WHERE terminal is proposed for Quack.
+
+## Open
+
+- Pack each block's B once per graph version and re-measure; this would
+  confirm or refute finding 4.
+- A locality-preserving G2L order (for example a bandwidth-reducing order) and
+  its effect on the block count.
+- A sparse-block realization (CSR-in-block or mask-fold) for low-fill blocks,
+  with the crossover chosen per block.
+- The blasgraph HDR-semiring adapter: a named seam, not wired.
+- Every number above misses the ms gate. The work is real dense block math,
+  but the gate is not met and this is not a pass.

.claude/board/entries/README.md

@@ -25,10 +25,11 @@ index row, (3) no duplicate entry id. Checks 1 and 2 are deliberately
 opposite directions; the stranding this convention prevents shows up in
 exactly one of them, never both.
 
-587 entries, 2026-08-06 .. 2026-11-08.
+588 entries, 2026-08-06 .. 2026-11-08.
 
 | date | entry id | finding | file |
 |---|---|---|---|
+| 2026-11-08 | `dr-virtual-csr-bf16-gemm` |  | [2026-11-08-dr-virtual-csr-bf16-gemm.md](2026-11-08-dr-virtual-csr-bf16-gemm.md) |
 | 2026-11-08 | `dr-isa-stream-lab` |  | [2026-11-08-dr-isa-stream-lab.md](2026-11-08-dr-isa-stream-lab.md) |
 | 2026-11-07 | `D-DT64-7-L5b` |  | [2026-11-07-dt64-l5b-historical-basis.md](2026-11-07-dt64-l5b-historical-basis.md) |
 | 2026-11-07 | `D-DT64-7-L5a` |  | [2026-11-07-dt64-l5a-exact-basis-witness.md](2026-11-07-dt64-l5a-exact-basis-witness.md) |

.github/workflows/dr-isa-stream-lab.yml

@@ -69,5 +69,5 @@ jobs:
         run: cargo clippy --manifest-path crates/dr-isa-stream-lab/Cargo.toml --all-targets -- -D warnings
       - name: tests (prints the GEMM tier)
         run: cargo test --release --manifest-path crates/dr-isa-stream-lab/Cargo.toml -- --nocapture
-      - name: stage costs (every timed result is checked against the oracle)
-        run: cargo run --release --manifest-path crates/dr-isa-stream-lab/Cargo.toml --bin bench
+      - name: benchmark of every arm (each result checked against its oracle)
+        run: cargo run --release --manifest-path crates/dr-isa-stream-lab/Cargo.toml --bin bench_all

crates/dr-isa-stream-lab/BENCH.md

@@ -0,0 +1,45 @@
+# dr-isa-stream-lab: benchmark of every arm
+
+Generated by `cargo run --release --bin bench_all`. Every result is checked against its oracle before timing and on every timed run.
+
+Host: Granite Rapids VM, 4 vCPU, rustc 1.99.0. lance-graph `7c301f5` + this branch; ndarray `f6c1a07` (AdaWorldAPI/ndarray#414). Single-threaded throughout. This is not certification on dedicated silicon: AMX tile state can corrupt under host contention, which is why every timed result is checked.
+
+The lance-graph gate is ns per operation and µs per batch. Rows in milliseconds fail it and are recorded as findings, not passes.
+
+Host tile tier: `AMX TDPBF16PS`. n = 65536, 16 lanes, 21 warm runs (median).
+
+| Workload | Arm | Time | Notes |
+|---|---|---|---|
+| WHERE 16×16 | bind (cold) | 993.0 ns | DT64 lane binding, 16 fields |
+| WHERE 16×16 | 16 WHERE counts (popcount) | 2163.0 ns | the cube's fold terminal |
+| WHERE 16×16 | stage selection B (cold) | 713.8 µs | once per mask set, 2 MiB |
+| WHERE 16×16 | tile GEMM, whole lane, A in place | 172.5 µs | 1 call, A staged 0 B |
+| WHERE 16×16 | tile GEMM, whole lane, A staged (slab) | 394.7 µs | A copied once |
+| WHERE 16×16 | tile GEMM, whole lane, A staged (scattered) | 389.0 µs | lanes in separate allocations |
+| WHERE 16×16 | ndarray cblas_gemm_bf16bf16f32, whole 64K | 3397.3 µs | 1 call; A is the slab (pitch = K); B row-major |
+| WHERE 16×16 | reference: tile GEMM per 32-identity block | 2222.5 µs | 2,048 calls, B restaged per block |
+| WHERE 16×16 | reference: cube per-bit sum fold | 13.5 ms | scalar gather |
+| graph deg 8, 1/16 far | build virtual CSR (cold) | 35.1 ms | nnz 522603, 5194 blocks (7.9% of tile pairs) |
+| graph deg 8, 1/16 far | tile GEMM, A in place | 32.3 ms | 24280 calls |
+| graph deg 8, 1/16 far | tile GEMM, A staged | 33.6 ms |  |
+| graph deg 8, 1/16 far | ndarray cblas_gemm_bf16bf16f32 | 580.9 ms | 5194 calls, A staged (no lda) |
+| graph deg 8, 1/16 far | scalar CSR over the same lanes | 11.6 ms | non-GEMM baseline |
+| graph deg 8, all far | build virtual CSR (cold) | 182.3 ms | nnz 522981, 41756 blocks (63.7% of tile pairs) |
+| graph deg 8, all far | tile GEMM, A in place | 105.9 ms | 82647 calls |
+| graph deg 8, all far | tile GEMM, A staged | 120.2 ms |  |
+| graph deg 8, all far | ndarray cblas_gemm_bf16bf16f32 | 4730.9 ms | 41756 calls, A staged (no lda) |
+| graph deg 8, all far | scalar CSR over the same lanes | 17.4 ms | non-GEMM baseline |
+| graph deg 64, own tile | build virtual CSR (cold) | 108.2 ms | nnz 4194304, 256 blocks (0.4% of tile pairs) |
+| graph deg 64, own tile | tile GEMM, A in place | 15.3 ms | 4096 calls |
+| graph deg 64, own tile | tile GEMM, A staged | 12.5 ms |  |
+| graph deg 64, own tile | ndarray cblas_gemm_bf16bf16f32 | 51.4 ms | 256 calls, A staged (no lda) |
+| graph deg 64, own tile | scalar CSR over the same lanes | 75.0 ms | non-GEMM baseline |
+
+## Reading
+
+- **The 256 × 256 cube is the whole universe.** One lane or one WHERE mask over 65,536 identities is a single 256 × 256 matrix, and the WHERE rows run over all of it in one call per 16 lanes. The graph's edge relation is identity × identity: 256 × 256 tile pairs, each a 256 × 256 block. The virtual CSR addresses that second level by `(tile_r, tile_c, slot_r, slot_c)`.
+- **WHERE (selection about 50% dense).** The strided tile GEMM with A in place is the fastest arm: 172 µs for 16.7M MACs. Copying A first doubles the cost.
+- **ndarray `cblas_gemm_bf16bf16f32` is about 20× slower than the strided tile path on the same problem**, although A is the slab itself at pitch = K. `amx_matmul::matmul_bf16_to_f32` copies both operands (`pack_contig`) and calls `bf16_tile_gemm_16x16` per 16×16 tile, which VNNI-packs B and allocates on every call. Before this entry can serve a stream it needs a prepacked or strided-B variant.
+- **Graph aggregation.** A block GEMM's cost is set by block area, not by edges. It loses to a scalar CSR loop at about 0.15% fill (degree 8) and wins about 5–6× at 25% fill (degree 64 within the tile).
+- **Long-range edges fill the block layout.** One long-range edge per identity grows the block count about 8×. A virtual CSR needs a locality-preserving G2L order.
+- **Every graph row is in milliseconds.** None passes the gate.