review: gate gridlake_field_tile to x86_64 + use canonical simd surface (#233, codex P2)

claude · claude · commit 368e01d026a2 · 2026-07-04T10:44:04.000Z
The example imported the x86_64-only hpc::amx_matmul / hpc::bf16_tile_gemm paths under required-features = ["std"], so a cross-target `cargo check --examples` would fail on non-x86 with unresolved imports. Fixes: - Wrap the whole probe in `#[cfg(target_arch = "x86_64")] mod amx { ... }` with a `#[cfg(not(target_arch = "x86_64"))]` fallback `main` that prints the x86-only requirement. Non-x86 builds now compile (and skip) cleanly. - Switch to the canonical `ndarray::simd::*` surface (per the W1a consumer contract — consumers import from `ndarray::simd`, not raw `hpc::*`): `amx_available`, `bf16_tile_gemm_16x16_amx` (the TDPBF16PS tile wrapper, aliased to the body's name), `bf16_tile_gemm_tier` — all arch-gated re-exports. Behaviour unchanged on x86: still runs the real AMX TDPBF16PS tile op (amx_available = true), 0.095% / 0.401% Frobenius, bit-exact on integer operands. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MLBnPuScZy6w9di2QEjsXM
diff --git a/examples/gridlake_field_tile.rs b/examples/gridlake_field_tile.rs
@@ -30,195 +30,204 @@
 //!
 //! Run: `cargo run --release --example gridlake_field_tile --features std`
 
-use ndarray::hpc::amx_matmul::amx_available;
-use ndarray::hpc::bf16_tile_gemm::{bf16_tile_gemm_16x16, bf16_tile_gemm_tier};
-use ndarray::hpc::quantized::{f32_to_bf16_rounded, BF16};
-
-const M: usize = 16; // tile edge — the AMX tile / gridlake quadrant
-const KPAD: usize = 32; // bf16_tile_gemm requires K a multiple of 32; pad 16→32 with zeros
-
-fn mix(mut z: u64) -> u64 {
-    z = z.wrapping_add(0x9E37_79B9_7F4A_7C15);
-    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
-    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
-    z ^ (z >> 31)
-}
-
-/// Pack a row-major `M×KPAD` (or `KPAD×M`) f32 matrix into the `&[u16]` bf16
-/// operand the tile GEMM consumes (RNE, matching `VCVTNEPS2BF16`).
-fn to_bf16(src: &[f32]) -> Vec<u16> {
-    let mut b = vec![BF16(0); src.len()];
-    f32_to_bf16_rounded(src, &mut b);
-    b.iter().map(|x| x.0).collect()
-}
+// The AMX / BF16-tile path is x86_64-only: `hpc::amx_matmul` and
+// `hpc::bf16_tile_gemm` are `#[cfg(target_arch = "x86_64")]`, and the canonical
+// `ndarray::simd::*` re-exports of the tile ladder are arch-gated to match. Wrap
+// the whole probe in an x86 module with a non-x86 fallback `main`, so a
+// cross-target `cargo check --examples` does not hit unresolved imports.
+#[cfg(target_arch = "x86_64")]
+mod amx {
+    // Canonical public surface (`ndarray::simd::*`, per the W1a consumer contract),
+    // not the raw `hpc::*` paths. `bf16_tile_gemm_16x16_amx` is the tile-dispatching
+    // TDPBF16PS wrapper (aliased to the name the body already uses).
+    use ndarray::hpc::quantized::{f32_to_bf16_rounded, BF16};
+    use ndarray::simd::{amx_available, bf16_tile_gemm_16x16_amx as bf16_tile_gemm_16x16, bf16_tile_gemm_tier};
+
+    const M: usize = 16; // tile edge — the AMX tile / gridlake quadrant
+    const KPAD: usize = 32; // bf16_tile_gemm requires K a multiple of 32; pad 16→32 with zeros
+
+    fn mix(mut z: u64) -> u64 {
+        z = z.wrapping_add(0x9E37_79B9_7F4A_7C15);
+        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+        z ^ (z >> 31)
+    }
 
-/// `C[16×16] = A[16×16] · B[16×16]` via the SHIPPED BF16 tile op. A/B are padded
-/// on K from 16→32 with zeros (the pad contributes nothing), so this is one
-/// `bf16_tile_gemm_16x16` — one `TDPBF16PS` tile op on the AMX tier.
-fn tile_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f32; M * M] {
-    let mut a_pad = [0.0f32; M * KPAD];
-    let mut b_pad = [0.0f32; KPAD * M];
-    for i in 0..M {
-        for k in 0..M {
-            a_pad[i * KPAD + k] = a[i * M + k]; // A[16×32], cols 16..31 = 0
-            b_pad[k * M + i] = b[k * M + i]; // B[32×16], rows 16..31 = 0 (already)
-        }
+    /// Pack a row-major `M×KPAD` (or `KPAD×M`) f32 matrix into the `&[u16]` bf16
+    /// operand the tile GEMM consumes (RNE, matching `VCVTNEPS2BF16`).
+    fn to_bf16(src: &[f32]) -> Vec<u16> {
+        let mut b = vec![BF16(0); src.len()];
+        f32_to_bf16_rounded(src, &mut b);
+        b.iter().map(|x| x.0).collect()
     }
-    let (ab, bb) = (to_bf16(&a_pad), to_bf16(&b_pad));
-    let mut c = vec![0.0f32; M * M];
-    bf16_tile_gemm_16x16(&ab, &bb, &mut c, KPAD);
-    let mut out = [0.0f32; M * M];
-    out.copy_from_slice(&c);
-    out
-}
 
-/// Direct `f64` reference GEMM `C = A·B` (16×16).
-fn direct_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f64; M * M] {
-    let mut c = [0.0f64; M * M];
-    for i in 0..M {
-        for j in 0..M {
-            let mut s = 0.0f64;
+    /// `C[16×16] = A[16×16] · B[16×16]` via the SHIPPED BF16 tile op. A/B are padded
+    /// on K from 16→32 with zeros (the pad contributes nothing), so this is one
+    /// `bf16_tile_gemm_16x16` — one `TDPBF16PS` tile op on the AMX tier.
+    fn tile_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f32; M * M] {
+        let mut a_pad = [0.0f32; M * KPAD];
+        let mut b_pad = [0.0f32; KPAD * M];
+        for i in 0..M {
             for k in 0..M {
-                s += a[i * M + k] as f64 * b[k * M + j] as f64;
+                a_pad[i * KPAD + k] = a[i * M + k]; // A[16×32], cols 16..31 = 0
+                b_pad[k * M + i] = b[k * M + i]; // B[32×16], rows 16..31 = 0 (already)
             }
-            c[i * M + j] = s;
         }
+        let (ab, bb) = (to_bf16(&a_pad), to_bf16(&b_pad));
+        let mut c = vec![0.0f32; M * M];
+        bf16_tile_gemm_16x16(&ab, &bb, &mut c, KPAD);
+        let mut out = [0.0f32; M * M];
+        out.copy_from_slice(&c);
+        out
     }
-    c
-}
 
-fn transpose(a: &[f32; M * M]) -> [f32; M * M] {
-    let mut t = [0.0f32; M * M];
-    for i in 0..M {
-        for j in 0..M {
-            t[j * M + i] = a[i * M + j];
+    /// Direct `f64` reference GEMM `C = A·B` (16×16).
+    fn direct_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f64; M * M] {
+        let mut c = [0.0f64; M * M];
+        for i in 0..M {
+            for j in 0..M {
+                let mut s = 0.0f64;
+                for k in 0..M {
+                    s += a[i * M + k] as f64 * b[k * M + j] as f64;
+                }
+                c[i * M + j] = s;
+            }
         }
+        c
     }
-    t
-}
 
-/// `(max_abs_err, frobenius_relative_err)` of a tile-GEMM result vs the f64
-/// direct. Frobenius-relative `‖tile − direct‖_F / ‖direct‖_F` is the honest
-/// aggregate precision — a per-element max-relative blows up on the near-zero
-/// entries that a cancellation-heavy product (e.g. `Mᵀ·Σ·M`) naturally has,
-/// which measures the cancellation, not the kernel.
-fn errors(tile: &[f32; M * M], direct: &[f64; M * M]) -> (f64, f64) {
-    let mut max_abs = 0.0f64;
-    let mut num = 0.0f64; // ‖E‖_F²
-    let mut den = 0.0f64; // ‖direct‖_F²
-    for i in 0..M * M {
-        let e = tile[i] as f64 - direct[i];
-        max_abs = max_abs.max(e.abs());
-        num += e * e;
-        den += direct[i] * direct[i];
+    fn transpose(a: &[f32; M * M]) -> [f32; M * M] {
+        let mut t = [0.0f32; M * M];
+        for i in 0..M {
+            for j in 0..M {
+                t[j * M + i] = a[i * M + j];
+            }
+        }
+        t
     }
-    (max_abs, (num / den.max(1e-12)).sqrt())
-}
 
-fn main() {
-    println!("Gridlake field interdependency = BF16 16×16 AMX tile op — measured");
-    println!("  tile M=16 (the AMX tile; 64×64 gridlake = 4×4 = 16 of these), K padded 16→32");
-    println!(
-        "  shipped kernel tier on THIS host: {}  (amx_available = {})\n",
-        bf16_tile_gemm_tier(),
-        amx_available()
-    );
-
-    // ── field tile X: a correlated gridlake patch (smooth → strong interdependency)
-    let mut x = [0.0f32; M * M];
-    for r in 0..M {
-        for c in 0..M {
-            x[r * M + c] = 4.0 * (r as f32 * 0.4).sin() * (c as f32 * 0.35).cos() + 0.5 * (r as f32 - c as f32);
+    /// `(max_abs_err, frobenius_relative_err)` of a tile-GEMM result vs the f64
+    /// direct. Frobenius-relative `‖tile − direct‖_F / ‖direct‖_F` is the honest
+    /// aggregate precision — a per-element max-relative blows up on the near-zero
+    /// entries that a cancellation-heavy product (e.g. `Mᵀ·Σ·M`) naturally has,
+    /// which measures the cancellation, not the kernel.
+    fn errors(tile: &[f32; M * M], direct: &[f64; M * M]) -> (f64, f64) {
+        let mut max_abs = 0.0f64;
+        let mut num = 0.0f64; // ‖E‖_F²
+        let mut den = 0.0f64; // ‖direct‖_F²
+        for i in 0..M * M {
+            let e = tile[i] as f64 - direct[i];
+            max_abs = max_abs.max(e.abs());
+            num += e * e;
+            den += direct[i] * direct[i];
         }
+        (max_abs, (num / den.max(1e-12)).sqrt())
     }
 
-    // ── (1) intercorrelation C = Xᵀ·X — the inter-cell coupling as a field op ──
-    let xt = transpose(&x);
-    let corr_tile = tile_matmul(&xt, &x);
-    let corr_direct = direct_matmul(&xt, &x);
-    let (ca, cr) = errors(&corr_tile, &corr_direct);
-    println!("  (1) intercorrelation  C = Xᵀ·X   (one tile op)");
-    println!("      max_abs_err={ca:.4}  frobenius_rel_err={:.3}%  (BF16-precision class)", cr * 100.0);
-
-    // ── (2) covariance sandwich Σ' = Mᵀ·Σ·M — the splat EWA projection shape ──
-    // Σ = a symmetric PSD covariance (Xᵀ·X); M = a coupling/projection matrix.
-    let sigma_f: [f32; M * M] = {
-        let mut s = [0.0f32; M * M];
-        for i in 0..M * M {
-            s[i] = corr_direct[i] as f32;
+    pub fn run() {
+        println!("Gridlake field interdependency = BF16 16×16 AMX tile op — measured");
+        println!("  tile M=16 (the AMX tile; 64×64 gridlake = 4×4 = 16 of these), K padded 16→32");
+        println!(
+            "  shipped kernel tier on THIS host: {}  (amx_available = {})\n",
+            bf16_tile_gemm_tier(),
+            amx_available()
+        );
+
+        // ── field tile X: a correlated gridlake patch (smooth → strong interdependency)
+        let mut x = [0.0f32; M * M];
+        for r in 0..M {
+            for c in 0..M {
+                x[r * M + c] = 4.0 * (r as f32 * 0.4).sin() * (c as f32 * 0.35).cos() + 0.5 * (r as f32 - c as f32);
+            }
         }
-        s
-    };
-    let mut m_proj = [0.0f32; M * M];
-    for r in 0..M {
-        for c in 0..M {
-            m_proj[r * M + c] = if r == c {
-                1.0
-            } else {
-                0.15 * (r as f32 - c as f32).signum()
-            };
+
+        // ── (1) intercorrelation C = Xᵀ·X — the inter-cell coupling as a field op ──
+        let xt = transpose(&x);
+        let corr_tile = tile_matmul(&xt, &x);
+        let corr_direct = direct_matmul(&xt, &x);
+        let (ca, cr) = errors(&corr_tile, &corr_direct);
+        println!("  (1) intercorrelation  C = Xᵀ·X   (one tile op)");
+        println!("      max_abs_err={ca:.4}  frobenius_rel_err={:.3}%  (BF16-precision class)", cr * 100.0);
+
+        // ── (2) covariance sandwich Σ' = Mᵀ·Σ·M — the splat EWA projection shape ──
+        // Σ = a symmetric PSD covariance (Xᵀ·X); M = a coupling/projection matrix.
+        let sigma_f: [f32; M * M] = {
+            let mut s = [0.0f32; M * M];
+            for i in 0..M * M {
+                s[i] = corr_direct[i] as f32;
+            }
+            s
+        };
+        let mut m_proj = [0.0f32; M * M];
+        for r in 0..M {
+            for c in 0..M {
+                m_proj[r * M + c] = if r == c {
+                    1.0
+                } else {
+                    0.15 * (r as f32 - c as f32).signum()
+                };
+            }
         }
-    }
-    // tile path: T = Σ·M, then Σ' = Mᵀ·T  (two tile ops — the sandwich)
-    let t_tile = tile_matmul(&sigma_f, &m_proj);
-    let mt = transpose(&m_proj);
-    let sand_tile = tile_matmul(&mt, &t_tile);
-    // direct reference: Mᵀ·Σ·M in f64
-    let t_direct = {
-        let mut td = [0.0f32; M * M];
-        let d = direct_matmul(&sigma_f, &m_proj);
+        // tile path: T = Σ·M, then Σ' = Mᵀ·T  (two tile ops — the sandwich)
+        let t_tile = tile_matmul(&sigma_f, &m_proj);
+        let mt = transpose(&m_proj);
+        let sand_tile = tile_matmul(&mt, &t_tile);
+        // direct reference: Mᵀ·Σ·M in f64
+        let t_direct = {
+            let mut td = [0.0f32; M * M];
+            let d = direct_matmul(&sigma_f, &m_proj);
+            for i in 0..M * M {
+                td[i] = d[i] as f32;
+            }
+            td
+        };
+        let sand_direct = direct_matmul(&mt, &t_direct);
+        let (sa, sr) = errors(&sand_tile, &sand_direct);
+        println!("  (2) covariance sandwich  Σ' = Mᵀ·Σ·M   (two tile ops = the splat EWA projection)");
+        println!(
+            "      max_abs_err={sa:.4}  frobenius_rel_err={:.3}%  (== hpc::splat3d::sandwich_x16 shape;",
+            sr * 100.0
+        );
+        println!("      per-entry max-rel is meaningless here — Mᵀ·Σ·M cancels to near-zero entries)");
+
+        // ── (3) bit-exactness on bf16-exact integer operands (module guarantee) ──
+        let mut ai = [0.0f32; M * M];
+        let mut bi = [0.0f32; M * M];
         for i in 0..M * M {
-            td[i] = d[i] as f32;
+            ai[i] = ((mix(i as u64) % 16) as f32) - 8.0; // small ints, bf16-exact
+            bi[i] = ((mix(i as u64 ^ 0x5555) % 16) as f32) - 8.0;
         }
-        td
-    };
-    let sand_direct = direct_matmul(&mt, &t_direct);
-    let (sa, sr) = errors(&sand_tile, &sand_direct);
-    println!("  (2) covariance sandwich  Σ' = Mᵀ·Σ·M   (two tile ops = the splat EWA projection)");
-    println!(
-        "      max_abs_err={sa:.4}  frobenius_rel_err={:.3}%  (== hpc::splat3d::sandwich_x16 shape;",
-        sr * 100.0
-    );
-    println!("      per-entry max-rel is meaningless here — Mᵀ·Σ·M cancels to near-zero entries)");
-
-    // ── (3) bit-exactness on bf16-exact integer operands (module guarantee) ──
-    let mut ai = [0.0f32; M * M];
-    let mut bi = [0.0f32; M * M];
-    for i in 0..M * M {
-        ai[i] = ((mix(i as u64) % 16) as f32) - 8.0; // small ints, bf16-exact
-        bi[i] = ((mix(i as u64 ^ 0x5555) % 16) as f32) - 8.0;
-    }
-    let int_tile = tile_matmul(&ai, &bi);
-    let int_direct = direct_matmul(&ai, &bi);
-    let bit_exact = (0..M * M).all(|i| int_tile[i] as f64 == int_direct[i]);
-    println!("  (3) bit-exact on bf16-exact integer operands (|Σ| < 2^24): {bit_exact}");
-
-    // ── (4) throughput of the 16×16 tile op on this host ──
-    let iters = 200_000usize;
-    let t0 = std::time::Instant::now();
-    let mut acc = 0.0f32;
-    for it in 0..iters {
-        // vary the operand slightly so nothing is optimized away
-        let mut xv = x;
-        xv[0] += (it & 0x7) as f32 * 0.01;
-        let c = tile_matmul(&xv, &m_proj);
-        acc += c[0] + c[M * M - 1];
-    }
-    let dt = t0.elapsed().as_secs_f64();
-    let tiles_s = iters as f64 / dt;
-    println!(
-        "  (4) throughput: {iters} tile ops in {:.1} ms → {:.2} M tile-ops/s  (checksum {acc:.1})",
-        dt * 1000.0,
-        tiles_s / 1e6
-    );
-
-    // hard gate: BF16-precision class (rel err small) + integer bit-exactness
-    assert!(cr < 0.05, "intercorrelation rel err too high for BF16 class");
-    assert!(sr < 0.05, "sandwich rel err too high for BF16 class");
-    assert!(bit_exact, "bf16-exact integer tile op is NOT bit-exact vs direct");
-
-    println!(
-        "\n  MEASURED CONCLUSION: YES — a gridlake field interdependency reduces to a\n\
+        let int_tile = tile_matmul(&ai, &bi);
+        let int_direct = direct_matmul(&ai, &bi);
+        let bit_exact = (0..M * M).all(|i| int_tile[i] as f64 == int_direct[i]);
+        println!("  (3) bit-exact on bf16-exact integer operands (|Σ| < 2^24): {bit_exact}");
+
+        // ── (4) throughput of the 16×16 tile op on this host ──
+        let iters = 200_000usize;
+        let t0 = std::time::Instant::now();
+        let mut acc = 0.0f32;
+        for it in 0..iters {
+            // vary the operand slightly so nothing is optimized away
+            let mut xv = x;
+            xv[0] += (it & 0x7) as f32 * 0.01;
+            let c = tile_matmul(&xv, &m_proj);
+            acc += c[0] + c[M * M - 1];
+        }
+        let dt = t0.elapsed().as_secs_f64();
+        let tiles_s = iters as f64 / dt;
+        println!(
+            "  (4) throughput: {iters} tile ops in {:.1} ms → {:.2} M tile-ops/s  (checksum {acc:.1})",
+            dt * 1000.0,
+            tiles_s / 1e6
+        );
+
+        // hard gate: BF16-precision class (rel err small) + integer bit-exactness
+        assert!(cr < 0.05, "intercorrelation rel err too high for BF16 class");
+        assert!(sr < 0.05, "sandwich rel err too high for BF16 class");
+        assert!(bit_exact, "bf16-exact integer tile op is NOT bit-exact vs direct");
+
+        println!(
+            "\n  MEASURED CONCLUSION: YES — a gridlake field interdependency reduces to a\n\
          \x20 BF16 16×16 tile op (one `bf16_tile_gemm_16x16`, i.e. one `TDPBF16PS` on\n\
          \x20 the AMX tier), matching the direct f64 reference to BF16 precision and\n\
          \x20 BIT-EXACT for bf16-exact integer operands. It IS the same primitive as:\n\
@@ -228,5 +237,16 @@ fn main() {
          \x20 Interdependency / intercorrelation / perturbation / transform / covariance\n\
          \x20 projection are ONE op: the 16×16 BF16 tile GEMM. The 64×64 gridlake is\n\
          \x20 4×4 = 16 of these tiles."
-    );
+        );
+    }
+} // mod amx
+
+#[cfg(target_arch = "x86_64")]
+fn main() {
+    amx::run();
+}
+
+#[cfg(not(target_arch = "x86_64"))]
+fn main() {
+    eprintln!("gridlake_field_tile requires x86_64 (AMX TDPBF16PS / AVX-512 VDPBF16PS BF16 tile GEMM).");
 }
