diff --git a/.claude/board/ISSUES.md b/.claude/board/ISSUES.md index e371fcf4e..1934981fc 100644 --- a/.claude/board/ISSUES.md +++ b/.claude/board/ISSUES.md @@ -1,3 +1,12 @@ +## ISS-BF16-POLYFILL-PER-ELEMENT-ALLOC — ndarray's BF16 16×16 polyfill allocates per output element (2026-10-20) + +**Status:** OPEN. **Basis:** VERIFIED-IN-CODE (`ndarray/src/simd_ops.rs`, `bf16_tile_gemm_16x16`) and MEASURED (D-RMAP-2 §5a: a 256 × 256 `A·A` takes about 19 ms, about 0.9 GMAC/s, on an AVX-512 host without AMX or `avx512_bf16`). + +- Per 16 × 16 call: two decode buffers (`vec!` for A and B), then for each of the 256 output elements a fresh `vec![0.0; k]` filled by a strided read of a B column, an FMA chain and a horizontal `reduce_sum`. About 257 heap allocations per call. +- `bf16_tile_gemm_16x16_packed` on the same tier unpacks the VNNI buffer back to row-major (another allocation) before calling it. +- Consequence: on hosts without the tile tiers, the BF16 tile path cannot be compared with the bit arm on cost; it measures allocation, not arithmetic. +- **What closes it:** an outer-product polyfill in ndarray (broadcast `A[i, kk]`, FMA into a 16-lane accumulator per row, no per-element buffer), with the existing exactness tests, then the D-RMAP-2 §5a probe re-run. + ## ISS-MORK-LIBRARY-GLOBAL-ALLOCATOR — linking MORK replaces the binary's allocator (2026-10-19) **Status:** OPEN. **Basis:** MEASURED (`crates/mork-pathmap-probe`, a `#[global_allocator]` in the probe fails to link: "conflicts with global allocator in: pathmap"). diff --git a/.claude/board/STATUS_BOARD.md b/.claude/board/STATUS_BOARD.md index 526bb675e..82db799dd 100644 --- a/.claude/board/STATUS_BOARD.md +++ b/.claude/board/STATUS_BOARD.md @@ -40,7 +40,7 @@ These PRs merged with no board rows. Every row below is a PROPOSAL (documents an | D-RPR-2 | premise revision with selective retraction | BLOCKED on D-TGA-3 (per-derived-cell provenance) | stale-premise refusal only, until then | | D-RPR-3 | MeTTa-facing parity of a chain's answers | BLOCKED on D-ASM-4 (not authorized) | — | | D-MPS-8 | native join over a versioned relationship table: R and S laid out by `y` slot, leapfrog intersection as one mask-risc `And`, per-side generations and reseal, per-thought filter planes (`mork-pathmap-schnellverschluss-v1.md` §7b; `crates/mork-pathmap-probe/src/join.rs`) | MEASURED, probe only; not admitted (D-RMAP-1 BIND membrane review) | 8 join tests + native arms in the differential suite vs the oracle; 7 guards disable-verified; cyclic bodies open | -| D-RMAP-2 | context-declaration reuse and GraphBLAS/BF16 tile probe (`d-rmap-2-context-tile-algebra-probe-v1.md`) | PROPOSAL; not run | arms 1–5 in the plan; 297× (D-GTM-0j) is not evidence | +| D-RMAP-2 | context-declaration reuse and GraphBLAS/BF16 tile probe (`d-rmap-2-context-tile-algebra-probe-v1.md`) | LAB ONLY, not admitted to production; PARTIAL: one cell MEASURED (§5a, polyfill tier only); tile tiers and GraphBLAS arm not run | arms 1–5 in the plan; 297× (D-GTM-0j) is not evidence | | D-SRAM-0 | inventory Soufflé AST/ast2ram/RAM passes against D-ASM/D-RA/D-TGA owners | Queued | evidence per mapped pass; no code imported | | D-SRAM-1 | bound AST → RAM-like schedule for pure conjunctive selection | Queued | same answers and read sets on hand-built fixtures; unsupported ops typed-refused | | D-SRAM-2 | cold optimizer rewrites (SCC/strata, magic sets, ordering, dedup) | Queued (research until D-ASM AST placement and D-RA-2) | optimized == unoptimized == oracle; rewrite off removes only the plan shape; disabling a semantic guard fails | diff --git a/.claude/board/entries/2026-10-20-d-rmap-2-tile-gemm-cell.md b/.claude/board/entries/2026-10-20-d-rmap-2-tile-gemm-cell.md new file mode 100644 index 000000000..c45e2a9eb --- /dev/null +++ b/.claude/board/entries/2026-10-20-d-rmap-2-tile-gemm-cell.md @@ -0,0 +1,13 @@ +# D-RMAP-2 tile arm: Boolean composition and triangles on one focus (2026-10-20) + +**Status:** MEASURED (one cell of D-RMAP-2), probe only, one host. Record: `.claude/plans/d-rmap-2-context-tile-algebra-probe-v1.md` §5a. Code: `crates/mork-pathmap-probe/src/tile_gemm.rs`. + +## Measured + +- On a 256 × 256 focus, `ndarray::simd::bf16_tile_gemm_16x16_packed` gives walk counts equal to a triple-loop oracle exactly, and `count > 0` equals Boolean composition (1 %, 10 %, 50 %, skewed hubs). 5 guards disable-verified. +- The tier that ran is the F32x16 FMA polyfill (no AMX, no `avx512_bf16` on this host): `A·A` about 19 ms. The bit arm: composition 2.6–25 µs, walk counts by popcount 165–192 µs, triangles 1.9–72 µs. Staging bits to BF16 (16 KiB to 256 KiB) costs 150–400 µs, more than the whole bit composition. +- The polyfill allocates per output element (`ISS-BF16-POLYFILL-PER-ELEMENT-ALLOC`), so this host measures allocation, not GEMM. + +## Open + +- AMX and `VDPBF16PS` hosts; blocked tiles above one focus; the GraphBLAS arm. diff --git a/.claude/board/entries/README.md b/.claude/board/entries/README.md index 19ed0f174..434f88182 100644 --- a/.claude/board/entries/README.md +++ b/.claude/board/entries/README.md @@ -25,13 +25,14 @@ 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. -377 entries, 2026-08-06 .. 2026-10-20. +378 entries, 2026-08-06 .. 2026-10-20. | date | entry id | finding | file | |---|---|---|---| | 2026-10-20 | `relstore-quack-petgraph-plan` | | [2026-10-20-relstore-quack-petgraph-plan.md](2026-10-20-relstore-quack-petgraph-plan.md) | | 2026-10-20 | `D-RST` | | [2026-10-20-relstore-geometry-not-graph-api.md](2026-10-20-relstore-geometry-not-graph-api.md) | | 2026-10-20 | `recursive-perspective-rebinding` | | [2026-10-20-recursive-perspective-rebinding.md](2026-10-20-recursive-perspective-rebinding.md) | +| 2026-10-20 | `D-RMAP-2` | | [2026-10-20-d-rmap-2-tile-gemm-cell.md](2026-10-20-d-rmap-2-tile-gemm-cell.md) | | 2026-10-19 | `post-merge-recovery-1601-1617` | | [2026-10-19-post-merge-recovery-1601-1617.md](2026-10-19-post-merge-recovery-1601-1617.md) | | 2026-10-19 | `native-join-relationship-table` | | [2026-10-19-native-join-relationship-table.md](2026-10-19-native-join-relationship-table.md) | | 2026-10-19 | `mork-pathmap-schnellverschluss` | | [2026-10-19-mork-pathmap-schnellverschluss.md](2026-10-19-mork-pathmap-schnellverschluss.md) | diff --git a/.claude/plans/d-rmap-2-context-tile-algebra-probe-v1.md b/.claude/plans/d-rmap-2-context-tile-algebra-probe-v1.md index ab4fe4721..e9ac17ca0 100644 --- a/.claude/plans/d-rmap-2-context-tile-algebra-probe-v1.md +++ b/.claude/plans/d-rmap-2-context-tile-algebra-probe-v1.md @@ -1,5 +1,7 @@ # D-RMAP-2: context declaration reuse and a GraphBLAS/BF16 tile probe +**LAB EXPLORATION ONLY (operator, 2026-10-20): nothing here is admitted to production, and no production code may be rewired to it. Unfit to promote until further testing.** + **Status: PROPOSAL / PREREGISTERED TEST DESIGN, 2026-10-19.** No runtime, ABI, dependency, opcode, CE64 interpretation or performance result is added. Extends [D-RMAP-1](d-rmap-implicit-relationship-register-v1.md) and its controlling [D-RMAP-1A ownership correction](d-rmap-separate-tables-thought-palettes-v1.md). @@ -92,6 +94,42 @@ Report: total bound cost, sorted-index creation/maintenance, native scan time, r - For sparse 1-hop, bitset or direct indexed seeks may defeat GEMM. For repeated dense 2-hop/all-pairs, tile GEMM may have a valid economic region. A loss in the full-cost test remains a loss. - Do not add an MKL/GraphBLAS C dependency, hot JSON representation, new runtime/AST, unauthorized gather or hidden cross-coordinate index hop. +## 5a. First measured cell: Boolean composition and triangles on one 256 × 256 focus (2026-10-20) + +> **Status: MEASURED, probe only, one host.** Code: `crates/mork-pathmap-probe/src/tile_gemm.rs`, `tests/tile_gemm.rs` (4 tests; 5 guards disable-verified), `examples/tile_gemm_cost.rs` (`cargo run --release --example tile_gemm_cost`). Release, warm, medians of 7, every timed answer checked against a triple-loop oracle. The GEMM arm is `ndarray::simd::bf16_tile_gemm_16x16_packed` through the facade, not a new kernel. + +**Tier that ran: the F32x16 FMA polyfill** (`bf16_tile_gemm_tier()`). This host has `avx512_vnni` and `avx512bw` but neither AMX nor `avx512_bf16`. No tile-tier number exists here; the AMX and `VDPBF16PS` cells of §4 remain unmeasured. This is the hardware, not detection (ndarray `.claude/AMX_GOTCHAS.md` Gotcha 4: a wrong `arch_prctl` syscall number once made `amx_available()` false on real AMX hosts). Checked directly on this host: CPUID.7.0 EDX bits 22/24/25 (AMX-BF16/TILE/INT8) = 0, CPUID.7.1 EAX bit 5 (AVX512-BF16) = 0, XCR0 = `0xe7` (XTILECFG/XTILEDATA bits 17/18 clear), and `arch_prctl(ARCH_REQ_XCOMP_PERM, XFEATURE_XTILEDATA)` on syscall 158 returns `-95` (EOPNOTSUPP). The AMX status in the gotcha file was measured on an Emerald Rapids VM, which this one is not. + +| Map | edges | triangles | bit compose | bit walk counts (popcount) | bit triangles | BF16 stage | BF16 `A·A` (packed / row-major) | +| --- | ---: | ---: | ---: | ---: | ---: | ---: | ---: | +| 1 % | 658 | 23 | 3.6 µs | 175 µs | 1.9 µs | 155 µs | 19.6 / 19.4 ms | +| 10 % | 6 529 | 16 617 | 4.5 µs | 190 µs | 13.5 µs | 256 µs | 19.4 / 22.3 ms | +| 50 % | 32 776 | 2 099 676 | 25.2 µs | 192 µs | 71.8 µs | 400 µs | 19.5 / 18.9 ms | +| 1 % + 4 hubs | 2 706 | 20 273 | 2.6 µs | 165 µs | 5.8 µs | 150 µs | 19.7 / 19.4 ms | + +- **§3's exactness claim holds as measured.** BF16 walk counts equal the oracle exactly (`assert_eq!`, not a tolerance), and `count > 0` equals Boolean composition, on all four maps. +- **On this host the bit arm wins every cell**, including walk counts. Popcount of `row_i(A) ∧ row_j(Aᵀ)` gives the same counts as the GEMM in about 170–190 µs, against about 19 ms. Staging alone (bits to dense BF16, 16 KiB to 256 KiB) costs more than the whole bit composition. +- **The 19 ms is the polyfill's shape, not GEMM's.** `ndarray::simd_ops::bf16_tile_gemm_16x16` allocates two decode buffers per call, then one strided B-column buffer and one horizontal reduce per output element: about 257 heap allocations per 16 × 16 call, 65 000 per 256 × 256 product. The packed entry on this tier also unpacks and allocates per call. Recorded as `ISS-BF16-POLYFILL-PER-ELEMENT-ALLOC`; the row-major entry is no faster, so the unpack is not the cost. +- **Two conclusions, kept apart.** (1) *Measured:* on this host's F32x16 polyfill tier, the bit arm wins every cell, and BF16 is exact. (2) *Not measured:* nothing here says anything about BF16 or AMX tile hardware performance. The polyfill's cost is dominated by allocation, so it is not even a lower bound on what `TDPBF16PS` or `VDPBF16PS` would cost. No layout decision may cite (1) as evidence about (2). +- **Not measured:** AMX and `VDPBF16PS` hosts; dense relations above one focus (blocked 256 × 256 tiles); weighted or evidence-bearing relations, which §3 excludes anyway. + +## 5b. Tile GEMM is a consumer of prepared relationships, not a gather solution (operator, 2026-10-20) + +> **Status: RULING-SHAPED CLARIFICATION, recorded as provenance.** Shared with D-RST (`relstore-quack-petgraph-traversal-v1.md` §5b), which carries the cost model and the falsifier F8. + +Tile GEMM speeds up arithmetic once relationship operands are in a compatible dense layout. It does not resolve global identities, align coordinate domains, scatter, gather or traverse sparse data. Those happen before the kernel and are paid there. A kernel loads its operands from memory, but only from tiles that something already built. + +**Three costs, always reported apart:** preparation (identity resolution, coordinate alignment, index traversal, gathering, packing, tile construction, sealing), execution (bit operations, SIMD, BF16/AMX GEMM on prepared operands) and maintenance (version checks, invalidation, incremental updates, repacking, rebuilds). A tile path wins only if `T_prepare + N·T_exec + M·T_maintain` falls below `N·T_traverse` for the measured reuse count `N` and invalidation count `M`. Break-even counts come from measured values, never from an assumed kernel speedup, and a warm kernel is never compared with a cold traversal. + +**What §5a measured, in those terms.** +- *Execution:* bit composition, popcount counts, triangles and the BF16 product, each on operands already in one 256 × 256 focus. +- *Part of preparation:* the "BF16 stage" column (bits to dense BF16), and it alone exceeds the whole bit composition. +- *Not measured:* the preparation that obtains the focus (identity resolution and slot assignment, measured separately in D-RMAP-1A §8), maintenance (no invalidation was exercised), and any reuse count. So §5a ranks execution paths on prepared operands; it says nothing yet about an end-to-end relational operation. + +**Boolean and arithmetic answers are different questions.** Reachability needs existence (`And`/`Or`); the BF16 product returns path multiplicities. §5a checked that `count > 0` equals Boolean composition. Using a count as existence is an explicit derivation step, never an implicit one. Exactness: the operands here are 0/1 and the accumulator is F32, so one product is exact up to 2^24. A count fed back as a BF16 operand (a further hop, or a weighted relation) is exact only up to 256, the BF16 significand bound. Such a path checks its counts against that bound and refuses beyond it. + +**No new production kernel is added for this hypothesis.** The probe uses `ndarray::simd::bf16_tile_gemm_16x16{,_packed}` through the facade and the existing mask-risc/ndarray bit algebra. The Boolean tile already exists as a mask-risc presence plane and as the 256 × 256 tile in `rmap.rs`. Sparse data is never forced into tiles in order to use GEMM. + ## 6. Recommendation Use **lance-context GenericStore/SchemaSpec only as a candidate cold durable declaration precedent**, **petgraph::visit as an optional read-only compatibility/oracle facade**, **GraphBLAS as the Boolean algebra oracle**, and **MKL BF16 GEMM as an explicitly falsifiable physical tile accelerator**. The map bank remains in independent typed tables; OGAR provides vocabulary; BIND selects and borrows; RelAlg describes; existing Quack/R2IL/Mask-RISC remains the execution family; CE64 is opaque. diff --git a/crates/mork-pathmap-probe/Cargo.toml b/crates/mork-pathmap-probe/Cargo.toml index 6e446813b..7498f61d3 100644 --- a/crates/mork-pathmap-probe/Cargo.toml +++ b/crates/mork-pathmap-probe/Cargo.toml @@ -34,6 +34,10 @@ pathmap = { git = "https://github.com/AdaWorldAPI/PathMap", rev = "ef02591" } # The native executor for the D-RMAP focus probe (`src/rmap.rs`): the prepared # relationship operands are run as its programs, never re-implemented here. lance-graph-mask-risc = { path = "../lance-graph-mask-risc" } +# The BF16 16x16 tile GEMM for the D-RMAP-2 tile arm (`src/tile_gemm.rs`), +# through the `ndarray::simd` facade. Same coordinate as mask-risc's, so the +# graph keeps one ndarray package. +ndarray = { git = "https://github.com/AdaWorldAPI/ndarray", branch = "master", default-features = false, features = ["std"] } # MORK's workspace names PathMap as `path = "../PathMap/"`, outside its own # repository, so a git dependency on MORK cannot resolve it. This patch points diff --git a/crates/mork-pathmap-probe/examples/tile_gemm_cost.rs b/crates/mork-pathmap-probe/examples/tile_gemm_cost.rs new file mode 100644 index 000000000..2aacb39b8 --- /dev/null +++ b/crates/mork-pathmap-probe/examples/tile_gemm_cost.rs @@ -0,0 +1,98 @@ +//! D-RMAP-2 tile arm: Boolean composition and triangle counts on one +//! 256 x 256 focus map, bit matrix vs BF16 16x16 tile GEMM. +//! +//! Run: `cd crates/mork-pathmap-probe && cargo run --release --example tile_gemm_cost`. +//! Release, warm caches, medians of `REPS`. Every timed answer is checked +//! against the triple-loop oracle first. + +use std::hint::black_box; +use std::time::{Duration, Instant}; + +use mork_pathmap_probe::tile_gemm::{ + Bf16Staged, BitMat, N, fixture, oracle_square, oracle_triangles, threshold, tier, + triangles_from_square, +}; + +const REPS: usize = 7; + +fn time(mut f: impl FnMut() -> T) -> (Duration, T) { + let mut ts = Vec::with_capacity(REPS); + let mut last = None; + for _ in 0..REPS { + let t = Instant::now(); + let v = black_box(f()); + ts.push(t.elapsed()); + last = Some(v); + } + ts.sort(); + (ts[REPS / 2], last.unwrap()) +} + +fn us(d: Duration) -> f64 { + d.as_secs_f64() * 1e6 +} + +fn main() { + println!("bf16 tier: {} | release, medians of {REPS}", tier()); + println!( + "case, edges, triangles, transpose_us, bit_compose_us, bit_counts_us, bit_triangles_us, bf16_stage_us, bf16_square_packed_us, bf16_square_rowmajor_us, threshold_us, bf16_triangles_us, bit_bytes, bf16_bytes" + ); + for (label, p, hubs, seed) in [ + ("1%", 10u32, 0usize, 1u64), + ("10%", 100, 0, 2), + ("50%", 500, 0, 3), + ("1% + 4 hubs", 10, 4, 4), + ] { + let a = fixture(p, hubs, seed); + let want_sq = oracle_square(&a); + let want_t = oracle_triangles(&a); + + let (t_tr, at) = time(|| black_box(&a).transpose()); + let (t_bc, c) = time(|| black_box(&a).compose(black_box(&at), black_box(&a))); + let mut want_bits = BitMat::empty(); + for i in 0..N { + for j in 0..N { + if want_sq[i * N + j] > 0 { + want_bits.set(i, j); + } + } + } + assert_eq!(c, want_bits); + let (t_bq, bq) = time(|| black_box(&a).square_counts(black_box(&at))); + assert_eq!(bq, want_sq); + let (t_bt, tri) = time(|| black_box(&a).triangles(black_box(&at))); + assert_eq!(tri, want_t); + + let (t_st, staged) = time(|| Bf16Staged::stage(black_box(&a))); + let mut sq = vec![0f32; N * N]; + let (t_sq, _) = time(|| { + black_box(&staged).square(black_box(&mut sq)); + }); + assert!(sq.iter().zip(&want_sq).all(|(g, &w)| *g == w as f32)); + let mut sq2 = vec![0f32; N * N]; + let (t_sq2, _) = time(|| { + black_box(&staged).square_rowmajor(black_box(&mut sq2)); + }); + assert_eq!(sq2, sq); + let (t_th, cb) = time(|| threshold(black_box(&sq))); + assert_eq!(cb, want_bits); + let (t_gt, tri2) = time(|| triangles_from_square(black_box(&a), black_box(&sq))); + assert_eq!(tri2, want_t); + + println!( + "{label}, {}, {want_t}, {:.1}, {:.1}, {:.1}, {:.1}, {:.1}, {:.1}, {:.1}, {:.1}, {:.1}, {}, {}", + a.edges(), + us(t_tr), + us(t_bc), + us(t_bq), + us(t_bt), + us(t_st), + us(t_sq), + us(t_sq2), + us(t_th), + us(t_gt), + 2 * N * N / 8, + staged.bytes() + ); + } +} diff --git a/crates/mork-pathmap-probe/src/lib.rs b/crates/mork-pathmap-probe/src/lib.rs index 14e593970..4101bd14f 100644 --- a/crates/mork-pathmap-probe/src/lib.rs +++ b/crates/mork-pathmap-probe/src/lib.rs @@ -26,6 +26,7 @@ pub mod boundary; pub mod join; pub mod reference; pub mod rmap; +pub mod tile_gemm; pub use boundary::{Backend, BoundIndex, Refusal, bind, capability}; diff --git a/crates/mork-pathmap-probe/src/tile_gemm.rs b/crates/mork-pathmap-probe/src/tile_gemm.rs new file mode 100644 index 000000000..63535f4bf --- /dev/null +++ b/crates/mork-pathmap-probe/src/tile_gemm.rs @@ -0,0 +1,315 @@ +//! D-RMAP-2 tile arm: Boolean relationship composition and triangle counting +//! over one 256 × 256 focus map, as a bit matrix and as a BF16 16 × 16 tile +//! GEMM. +//! +//! STATUS: PROBE. Results only. +//! +//! - **Bit arm.** `A` is 256 rows of 256 bits (8 KiB). Composition +//! `C = A ∘ B` is the outer-product form: for every coordinate `k`, each row +//! `i` whose bit `k` is set ORs in row `k` of `B`. Every row is addressed by +//! its own coordinate (`i`, `k`); no position is taken from a stored value. +//! - **GEMM arm.** `A` staged as dense BF16 (0 or 1.0, 128 KiB), multiplied by +//! `ndarray::simd::bf16_tile_gemm_16x16_packed` in 16 × 16 output blocks with +//! `K = 256`, then thresholded. For 0/1 operands the products are exact, and +//! every count is at most 256, far below 2^24, so `count > 0` equals Boolean +//! composition (D-RMAP-2 §3). The counts themselves are walk counts, which +//! the bit arm cannot give without popcounts. +//! +//! Triangles: `T = Σ_{i,j,k} A[i,j]·A[j,k]·A[k,i]`, the directed 3-cycles +//! counted once per rotation. Both arms compute it, against a triple-loop +//! oracle. Nothing here applies to weighted, negated or evidence-bearing +//! relations; that is the boundary D-RMAP-2 §3 states. + +use ndarray::simd::{ + PackedBf16B, bf16_tile_gemm_16x16, bf16_tile_gemm_16x16_packed, bf16_tile_gemm_tier, +}; + +/// Focus side. +pub const N: usize = 256; +const W: usize = N / 64; +const BF16_ONE: u16 = 0x3F80; + +/// A Boolean relation over one focus: row `i`, bit `j` set when `A[i, j]`. +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct BitMat { + rows: Vec<[u64; W]>, +} + +impl BitMat { + pub fn empty() -> Self { + BitMat { + rows: vec![[0; W]; N], + } + } + + pub fn from_edges(edges: impl IntoIterator) -> Self { + let mut m = Self::empty(); + for (i, j) in edges { + m.set(i as usize, j as usize); + } + m + } + + pub fn set(&mut self, i: usize, j: usize) { + self.rows[i][j >> 6] |= 1 << (j & 63); + } + + pub fn get(&self, i: usize, j: usize) -> bool { + self.rows[i][j >> 6] >> (j & 63) & 1 == 1 + } + + pub fn row(&self, i: usize) -> &[u64; W] { + &self.rows[i] + } + + pub fn edges(&self) -> usize { + self.rows + .iter() + .flatten() + .map(|w| w.count_ones() as usize) + .sum() + } + + pub fn transpose(&self) -> BitMat { + let mut t = Self::empty(); + for i in 0..N { + for j in 0..N { + if self.get(i, j) { + t.set(j, i); + } + } + } + t + } + + /// `C = A ∘ B` (Boolean), outer-product form over the transpose of `A`: + /// for each `k`, the rows `i` with `A[i, k]` OR in `B[k]`. + pub fn compose(&self, at: &BitMat, b: &BitMat) -> BitMat { + let mut c = Self::empty(); + for k in 0..N { + let bk = b.rows[k]; + for (wi, &word) in at.rows[k].iter().enumerate() { + let mut w = word; + while w != 0 { + let i = wi * 64 + w.trailing_zeros() as usize; + let ci = &mut c.rows[i]; + for x in 0..W { + ci[x] |= bk[x]; + } + w &= w - 1; + } + } + } + c + } + + /// Walk counts `(A · B)[i, j] = popcount(row_i(A) ∧ row_j(Bᵀ))`, the bit + /// arm's answer to the counts GEMM gives. Row-major 256 × 256. + pub fn square_counts(&self, bt: &BitMat) -> Vec { + let mut c = vec![0u32; N * N]; + for i in 0..N { + let ri = self.rows[i]; + for j in 0..N { + let cj = bt.rows[j]; + let mut n = 0u32; + for x in 0..W { + n += (ri[x] & cj[x]).count_ones(); + } + c[i * N + j] = n; + } + } + c + } + + /// Directed 3-cycles: for every edge `(i, j)`, the `k` with `A[j, k]` and + /// `A[k, i]` is `row_j(A) ∧ row_i(Aᵀ)`; popcount and sum. + pub fn triangles(&self, at: &BitMat) -> u64 { + let mut t = 0u64; + for i in 0..N { + let ti = at.rows[i]; + for (wj, &word) in self.rows[i].iter().enumerate() { + let mut w = word; + while w != 0 { + let j = wj * 64 + w.trailing_zeros() as usize; + let rj = self.rows[j]; + for x in 0..W { + t += u64::from((rj[x] & ti[x]).count_ones()); + } + w &= w - 1; + } + } + } + t + } +} + +/// `A` staged for the GEMM: row-major BF16, plus each 16-column block of `A` +/// as a packed `B` operand (`K = 256`). +pub struct Bf16Staged { + a: Vec, + b_blocks: Vec, + /// The same blocks row-major, for the polyfill entry. On a host without + /// AMX or `avx512bf16`, the packed entry unpacks (and allocates) per call + /// before taking the polyfill, so the row-major entry is that tier's + /// honest cost. + b_rows: Vec>, +} + +impl Bf16Staged { + pub fn stage(m: &BitMat) -> Self { + let mut a = vec![0u16; N * N]; + for i in 0..N { + for j in 0..N { + if m.get(i, j) { + a[i * N + j] = BF16_ONE; + } + } + } + let b_rows: Vec> = (0..N / 16) + .map(|bj| { + let mut blk = vec![0u16; N * 16]; + for k in 0..N { + blk[k * 16..k * 16 + 16] + .copy_from_slice(&a[k * N + bj * 16..k * N + bj * 16 + 16]); + } + blk + }) + .collect(); + let b_blocks = b_rows.iter().map(|b| PackedBf16B::pack(b, N)).collect(); + Bf16Staged { + a, + b_blocks, + b_rows, + } + } + + /// Bytes held: the BF16 matrix and the packed blocks. + pub fn bytes(&self) -> usize { + self.a.len() * 2 + + self + .b_blocks + .iter() + .map(|b| b.data().len() * 2) + .sum::() + } + + /// Walk counts through the row-major polyfill entry (no pack, no unpack). + pub fn square_rowmajor(&self, c: &mut [f32]) { + assert_eq!(c.len(), N * N); + let mut blk = [0f32; 256]; + for bi in 0..N / 16 { + let a_rows = &self.a[bi * 16 * N..(bi + 1) * 16 * N]; + for (bj, b) in self.b_rows.iter().enumerate() { + blk.fill(0.0); + bf16_tile_gemm_16x16(a_rows, b, &mut blk, N); + for r in 0..16 { + let at = (bi * 16 + r) * N + bj * 16; + c[at..at + 16].copy_from_slice(&blk[r * 16..r * 16 + 16]); + } + } + } + } + + /// Walk counts `C = A · A` in f32, row-major 256 × 256, through the + /// tier-dispatched packed entry. + pub fn square(&self, c: &mut [f32]) { + assert_eq!(c.len(), N * N); + let mut blk = [0f32; 256]; + for bi in 0..N / 16 { + let a_rows = &self.a[bi * 16 * N..(bi + 1) * 16 * N]; + for (bj, b) in self.b_blocks.iter().enumerate() { + blk.fill(0.0); + bf16_tile_gemm_16x16_packed(a_rows, b, &mut blk); + for r in 0..16 { + let at = (bi * 16 + r) * N + bj * 16; + c[at..at + 16].copy_from_slice(&blk[r * 16..r * 16 + 16]); + } + } + } + } +} + +/// The tier ndarray will take on this host, for the report. +pub fn tier() -> &'static str { + bf16_tile_gemm_tier() +} + +/// Boolean composition from walk counts: `C[i, j] = counts[i, j] > 0`. +pub fn threshold(counts: &[f32]) -> BitMat { + let mut m = BitMat::empty(); + for i in 0..N { + for j in 0..N { + if counts[i * N + j] > 0.0 { + m.set(i, j); + } + } + } + m +} + +/// Triangles from walk counts: `Σ_{i,j} (A²)[j, i] · A[i, j]`. +pub fn triangles_from_square(m: &BitMat, sq: &[f32]) -> u64 { + let mut t = 0u64; + for i in 0..N { + for j in 0..N { + if m.get(i, j) { + t += sq[j * N + i] as u64; + } + } + } + t +} + +/// Independent oracle: a triple loop over the bits. +pub fn oracle_square(m: &BitMat) -> Vec { + let mut c = vec![0u32; N * N]; + for i in 0..N { + for k in 0..N { + if m.get(i, k) { + for j in 0..N { + if m.get(k, j) { + c[i * N + j] += 1; + } + } + } + } + } + c +} + +pub fn oracle_triangles(m: &BitMat) -> u64 { + let sq = oracle_square(m); + let mut t = 0u64; + for i in 0..N { + for j in 0..N { + if m.get(i, j) { + t += u64::from(sq[j * N + i]); + } + } + } + t +} + +/// A deterministic fixture: each edge with probability `p_permille / 1000`; +/// `hubs` rows (and columns) are dense, which gives a skewed degree +/// distribution. +pub fn fixture(p_permille: u32, hubs: usize, seed: u64) -> BitMat { + let mut s = seed; + let mut next = || { + s = s.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = s; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + }; + let mut m = BitMat::empty(); + for i in 0..N { + for j in 0..N { + let hub = i < hubs || j < hubs; + if hub || (next() % 1000) < u64::from(p_permille) { + m.set(i, j); + } + } + } + m +} diff --git a/crates/mork-pathmap-probe/tests/tile_gemm.rs b/crates/mork-pathmap-probe/tests/tile_gemm.rs new file mode 100644 index 000000000..6e874d3ad --- /dev/null +++ b/crates/mork-pathmap-probe/tests/tile_gemm.rs @@ -0,0 +1,87 @@ +//! The D-RMAP-2 tile arm against a triple-loop oracle. + +use mork_pathmap_probe::tile_gemm::{ + Bf16Staged, BitMat, N, fixture, oracle_square, oracle_triangles, threshold, + triangles_from_square, +}; + +const CASES: [(u32, usize, u64); 4] = [(10, 0, 1), (100, 0, 2), (500, 0, 3), (10, 4, 4)]; + +#[test] +fn bit_and_bf16_composition_equal_the_oracle() { + for (p, hubs, seed) in CASES { + let a = fixture(p, hubs, seed); + let at = a.transpose(); + let want = oracle_square(&a); + let mut want_bits = BitMat::empty(); + for i in 0..N { + for j in 0..N { + if want[i * N + j] > 0 { + want_bits.set(i, j); + } + } + } + assert_eq!(a.compose(&at, &a), want_bits, "bit arm p={p} hubs={hubs}"); + assert_eq!(a.square_counts(&at), want, "bit counts p={p} hubs={hubs}"); + let staged = Bf16Staged::stage(&a); + let mut sq = vec![0f32; N * N]; + staged.square(&mut sq); + // Exact walk counts, not only their support. + for (got, &w) in sq.iter().zip(&want) { + assert_eq!(*got, w as f32, "bf16 count p={p} hubs={hubs}"); + } + assert_eq!( + threshold(&sq), + want_bits, + "bf16 threshold p={p} hubs={hubs}" + ); + // Not degenerate: two hops reach something one hop does not. + assert_ne!(want_bits, a, "p={p}"); + } +} + +#[test] +fn triangle_counts_agree_across_arms() { + for (p, hubs, seed) in CASES { + let a = fixture(p, hubs, seed); + let at = a.transpose(); + let want = oracle_triangles(&a); + assert_eq!(a.triangles(&at), want, "bit arm p={p} hubs={hubs}"); + let mut sq = vec![0f32; N * N]; + Bf16Staged::stage(&a).square(&mut sq); + assert_eq!( + triangles_from_square(&a, &sq), + want, + "bf16 arm p={p} hubs={hubs}" + ); + } + // Silence twin: a DAG (edges only i < j) has no directed cycle. + let dag = BitMat::from_edges( + (0..N).flat_map(|i| (i + 1..(i + 4).min(N)).map(move |j| (i as u8, j as u8))), + ); + assert!(dag.edges() > 0); + assert_eq!(dag.triangles(&dag.transpose()), 0); + assert_eq!(oracle_triangles(&dag), 0); +} + +/// The counts carry multiplicity the Boolean result drops: some pair is +/// reached by more than one walk, so a threshold is a real loss of +/// information, not a no-op. +#[test] +fn walk_counts_exceed_one_where_support_does_not_show_it() { + let a = fixture(100, 0, 2); + let mut sq = vec![0f32; N * N]; + Bf16Staged::stage(&a).square(&mut sq); + assert!(sq.iter().any(|&c| c > 1.0)); +} + +/// An asymmetric relation, so a composition that forgets the transpose is +/// caught. +#[test] +fn composition_respects_direction() { + let a = BitMat::from_edges([(0, 1), (1, 2)]); + let c = a.compose(&a.transpose(), &a); + assert!(c.get(0, 2)); + assert!(!c.get(2, 0)); + assert_eq!(c.edges(), 1); +}