lkdvos · Oct 7, 2026 · Sep 30, 2026 · Sep 30, 2026 · Sep 30, 2026 · Sep 30, 2026
diff --git a/Project.toml b/Project.toml
@@ -9,28 +9,32 @@ SIMD = "fdea26ae-647d-5447-a871-4b548cad5224"
 StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
+VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8"
 
 [compat]
 Aqua = "0.8"
 Bumper = "0.7"
 InteractiveUtils = "1"
 LinearAlgebra = "1"
+ParallelTestRunner = "2.9"
 Random = "1"
 SIMD = "3"
 Strided = "2.6"
 StridedViews = "0.5"
 TensorOperations = "5.8"
 TupleTools = "1.6"
+VectorInterface = "0.4.1, 0.5, 0.6, 0.7"
 Test = "1"
 julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Bumper = "8ce10254-0962-460f-a3d8-1f77fea1446e"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "Strided", "Aqua", "Bumper", "InteractiveUtils"]
+test = ["Test", "Random", "Strided", "Aqua", "Bumper", "InteractiveUtils", "ParallelTestRunner"]
diff --git a/README.md b/README.md
@@ -64,8 +64,8 @@ plan = plan_contract(StridedView(C), StridedView(A), (1, 2, 3),
 execute!(plan, 1.0, 0.0)
 ```
 
-`plan_contract` also takes `workspace=`/`allocator=` keywords (e.g. for
-Bumper-backed buffers); see its docstring.
+`plan_contract` also takes an `allocator=` keyword (e.g. for Bumper-backed
+buffers); see its docstring.
 
 ## API
 
@@ -95,7 +95,7 @@ and `B` into panels and drives a microkernel over register tiles of `C`.
 | `src/packing/` | packed-panel formats and the `pack_a!`/`pack_b!` packers with contiguous fast paths |
 | `src/microkernels/` | the kernel interface and the scalar, SIMD, planar, 1m and fmaddsub kernels |
 | `src/planning/` | label classification, conjugation, kernel selection, cache blocking, `plan_contract` |
-| `src/execution/` | the workspace, the five-loop nest, the specialised paths and the tile-by-tile oracle |
+| `src/execution/` | the workspace, the barriers, the five-loop nest and the specialised paths |
 | `src/integrations/` | `QuasiStridedBackend` for TensorOperations.jl |
 
 `test/` mirrors this layout.
@@ -127,8 +127,8 @@ Implemented:
   `eltype(C)` once. `conjA`/`conjB` and each operand's `StridedView.op` are
   applied during packing.
 - Zero steady-state allocation on Julia >= 1.11 for `execute!` on a reused
-  plan and for `tensorcontract!` through the backend (which pools workspaces
-  per task). On Julia 1.10 `SIMDKernel`'s accumulator is not kept in
+  plan; `tensorcontract!` through the backend plans each call, so it allocates
+  the plan's workspace through the given allocator. On Julia 1.10 `SIMDKernel`'s accumulator is not kept in
   registers, so calls allocate; results are unaffected.
 
 Not implemented:

diff --git a/benchmark/README.md b/benchmark/README.md
@@ -0,0 +1,29 @@
+# Benchmarks
+
+Single-core scripts; each script's header has its options and outputs. All run
+from the repo root with `julia --project=benchmark benchmark/<script> [args]`,
+and every `bench_*.jl` takes `--smoke` (one rep on a few cases) to check that it runs.
+
+| Script | Use it to |
+|:-|:-|
+| `bench_to_suite.jl` | compare QuasiStridedBackend with StridedBLAS on the upstream TensorOperations suite (`:contract`, `:network`); the overall picture after any change |
+| `plot_bench_to_suite.jl` | plot a `bench_to_suite.csv`: per dtype, both backends' GFLOP/s and their time ratio against arithmetic intensity, with `--dtypes`/`--categories`/`--tags` filters |
+| `bench_kernels.jl` | measure every register tile of the kernel menu against the default choice and OpenBLAS, after kernel or kernel-selection changes |
+| `bench_blocking_model.jl` | check the cache-blocking model against a sweep of `(m_block, k_block, n_block)`, after blocking changes |
+| `bench_mixed.jl` | compare the mixed real/complex kernels with promotion and the real GEMM of equal FMA count |
+| `bench_degenerate.jl` | compare the dot, outer-product and unpacked-B paths with the five-loop nest and BLAS |
+| `probes/probe_ttfx.jl` | time to first call (load, first plans, `@tensor`, `contract!`); one fresh process per run |
+| `probes/probe_call_floor.jl` | per-call time and bytes on tiny problems (`tensorcontract!`, `contract!`, reused-plan `execute!`) |
+| `probes/probe_specialisations.jl` | count compiled specialisations of the hot functions after a fixed workload |
+
+On the cluster, one script per job, one job per microarchitecture:
+
+    sbatch --reservation=rocky8 --constraint=<icelake|genoa|rome> benchmark/submit.sh <script> [args]
+
+A same-node A/B of two revisions, alternating them over `$QS_AB_ROUNDS` rounds:
+
+    sbatch --reservation=rocky8 --constraint=<icelake|genoa|rome> benchmark/submit_ab.sh <rev_a> <rev_b> <script> [args]
+
+`benchmark/submit_bench_to_suite.sh` runs the full suite sweep plus its plots.
+Results land in `benchmark/results/` (gitignored). Timings on a shared
+workstation are indicative only.
diff --git a/benchmark/bench_blocking_model.jl b/benchmark/bench_blocking_model.jl
@@ -1,24 +1,25 @@
-# (mc, kc, nc) plateau sweep: the analytical blocking model and the fallback
-# row against a log-spaced grid, at the engine's default kernel.
+# (m_block, k_block, n_block) plateau sweep: the analytical blocking model and
+# the fallback row against a log-spaced grid, at the engine's default kernel.
 #
 #   julia --project=benchmark benchmark/bench_blocking_model.jl [--reps 21] [--smoke] [--outdir DIR]
 #
-# Per real dtype: the named rows (`model` = `_modelled_blocking`, `fallback` =
-# `_fallback_blocking`) and the `wide` grid on every harness shape, plus 1-D
-# `nc` and `mc` slices through the model row on a 2048x256x2048 shape (the only
-# one large enough for `nc` to bind). Per complex dtype: the named rows scaled
-# by `_scale_blocking`. Every configuration of a shape is timed back to back.
+# Per real dtype: the named rows (`model` = `default_blocking` on this host,
+# `fallback` = on an undetected one) and the `wide` grid on every harness
+# shape, plus 1-D `n_block` and `m_block` slices through the model row on a
+# 2048x256x2048 shape (the only one large enough for `n_block` to bind). Per
+# complex dtype: the named rows. Every configuration of a shape is timed back to back.
 # The score of a configuration is the geomean over shapes of its time over the
 # best time at that shape. `--smoke` runs 1 rep on a few points.
 
 include(joinpath(@__DIR__, "harness.jl"))
 
-using QuasiStrided: Blocking, target_profile, _default_kernel, _modelled_blocking,
-    _fallback_blocking, _scale_blocking, complex_method
+using QuasiStrided: Blocking, target_profile, unknown_target, default_blocking,
+    kernel_from_shape, derived_shape
 
-const SMOKE = hasflag("smoke")
-const REPS = SMOKE ? 1 : argopt("reps", 21)
-thin(v) = SMOKE ? v[1:min(end, 3)] : v
+const REPS = reps_arg(21)
+const DTYPES = (Float64, Float32)
+const CDTYPES = (ComplexF64, ComplexF32)
+full_grid(m_blocks, k_blocks, n_blocks) = [(m_block, k_block, n_block) for k_block in k_blocks for m_block in m_blocks for n_block in n_blocks]
 
 const GRID_SHAPES = vcat(MAIN_SHAPES, EXTRA_SHAPES, SMALL_SHAPES)
 const BIG_SHAPE = ShapeSpec(2048, 256, 2048)
@@ -34,27 +35,29 @@ const OUTDIR = outdir()
 const CSV_PATH = joinpath(OUTDIR, "blocking_model.csv")
 const SUMMARY_PATH = joinpath(OUTDIR, "summary_blocking_model.txt")
 
-tup(b::Blocking) = (b.mc, b.kc, b.nc)
+tup(b::Blocking) = (b.m_block, b.k_block, b.n_block)
 named_points(rows) = [(string(k), tup(v)) for (k, v) in pairs(rows) if v !== nothing]
-named_rows(::Type{T}) where {T <: Real} = (model = _modelled_blocking(PROFILE, T), fallback = _fallback_blocking(T))
+detected = PROFILE.l1d.bytes > 0 && PROFILE.l2.bytes > 0
+named_rows(kernel) = (model = detected ? default_blocking(kernel, PROFILE) : nothing, fallback = default_blocking(kernel, unknown_target()))
+host_kernel(::Type{T}) where {T} = kernel_from_shape(derived_shape(PROFILE, T), T)
 
 csv = open(CSV_PATH, "w")
-println(csv, "set,dtype,shape,Ma,Ka,Na,mc,kc,nc,mc_eff,kc_eff,nc_eff,reps,median_seconds,gflops")
+println(csv, "set,dtype,shape,Ma,Ka,Na,m_block_requested,k_block_requested,n_block_requested,m_block,k_block,n_block,reps,median_seconds,gflops")
 
 function sweep_shape!(raw, kernel, ::Type{T}, spec, points) where {T}
     fx = build_plain(T, spec, MersenneTwister(0xB10C))
-    for (set, (mc, kc, nc)) in points
+    for (set, (m_block, k_block, n_block)) in points
         plan = plan_contract(
-            fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC; kernel = kernel, mc = mc, kc = kc, nc = nc
+            fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC; kernel = kernel, m_block = m_block, k_block = k_block, n_block = n_block
         )
         t = median_time_s(() -> execute!(plan, one(T), zero(T)); reps = REPS)
         b = plan.blocking
         println(
-            csv, "$set,$T,$(spec.name),$(spec.Ma),$(spec.Ka),$(spec.Na),$mc,$kc,$nc,$(b.mc),$(b.kc),$(b.nc),$REPS,",
+            csv, "$set,$T,$(spec.name),$(spec.Ma),$(spec.Ka),$(spec.Na),$m_block,$k_block,$n_block,$(b.m_block),$(b.k_block),$(b.n_block),$REPS,",
             @sprintf("%.9f,%.4f", t, gflops(T, spec.Ma, spec.Ka, spec.Na, t))
         )
         flush(csv)
-        push!(raw, (; set, dtype = T, shape = spec.name, mc, kc, nc, t))
+        push!(raw, (; set, dtype = T, shape = spec.name, m_block, k_block, n_block, t))
     end
     return nothing
 end
@@ -80,8 +83,8 @@ function summarize(io, raw, canaries)
     @printf(io, "canary spread = %.1f%%  %s\n", 100 * relative_spread(canaries), canaries)
     grid_names = [s.name for s in GRID_SHAPES]
     for T in DTYPES
-        println(io, "\n## $T   ", named_rows(T))
-        w = score(raw, T, grid_names, r -> (r.mc, r.kc, r.nc); sets = ("wide",))
+        println(io, "\n## $T   ", named_rows(host_kernel(T)))
+        w = score(raw, T, grid_names, r -> (r.m_block, r.k_block, r.n_block); sets = ("wide",))
         @printf(
             io, "wide grid: best %.4f %s  worst %.4f %s  spread %.1f%%\n",
             w[1][2], w[1][1], w[end][2], w[end][1], 100 * (w[end][2] / w[1][2] - 1)
@@ -92,7 +95,7 @@ function summarize(io, raw, canaries)
         for (set, g) in score(raw, T, grid_names, r -> r.set; sets = ("model", "fallback"))
             @printf(io, "  %-10s %.4f\n", set, g)
         end
-        for (s, f) in (("nslice", r -> r.nc), ("mslice", r -> r.mc))
+        for (s, f) in (("nslice", r -> r.n_block), ("mslice", r -> r.m_block))
             pts = sort([r for r in raw if r.dtype == T && r.set == s]; by = f)
             isempty(pts) && continue
             tb = minimum(r.t for r in pts)
@@ -114,26 +117,26 @@ raw = NamedTuple[]
 crng = MersenneTwister(0x0CA9A121)
 canaries = [run_canary(crng, "start")]
 for T in DTYPES
-    kernel = _default_kernel(T)
-    rows = named_rows(T)
-    println("\n$T kernel $(mr(kernel))x$(nr(kernel))/W$(lanewidth(kernel))  ", rows)
+    kernel = host_kernel(T)
+    rows = named_rows(kernel)
+    println("\n$T kernel $(tile_size(kernel, 1))x$(tile_size(kernel, 2))/W$(lanewidth(kernel))  ", rows)
     rows.model === nothing && @warn "no cache geometry detected: model undefined"
     grid = [("wide", c) for c in thin(WIDE[T])]
     for spec in thin(GRID_SHAPES)
         sweep_shape!(raw, kernel, T, spec, vcat(named_points(rows), grid))
     end
     m = something(rows.model, rows.fallback)
     slices = vcat(
-        named_points(rows), [("nslice", (m.mc, m.kc, n)) for n in thin(NSLICE)],
-        [("mslice", (x, m.kc, m.nc)) for x in thin(MSLICE)],
+        named_points(rows), [("nslice", (m.m_block, m.k_block, n)) for n in thin(NSLICE)],
+        [("mslice", (x, m.k_block, m.n_block)) for x in thin(MSLICE)],
     )
     sweep_shape!(raw, kernel, T, BIG_SHAPE, slices)
     push!(canaries, run_canary(crng, "after-$T"))
 end
 for T in CDTYPES
-    kernel = _default_kernel(T)
-    rows = map(v -> v === nothing ? nothing : _scale_blocking(v, complex_method(kernel)), named_rows(real(T)))
-    println("\n$T kernel $(mr(kernel))x$(nr(kernel))/W$(lanewidth(kernel))  ", rows)
+    kernel = host_kernel(T)
+    rows = named_rows(kernel)
+    println("\n$T kernel $(tile_size(kernel, 1))x$(tile_size(kernel, 2))/W$(lanewidth(kernel))  ", rows)
     for spec in thin(MAIN_SHAPES)
         sweep_shape!(raw, kernel, T, spec, named_points(rows))
     end

diff --git a/benchmark/bench_degenerate.jl b/benchmark/bench_degenerate.jl
@@ -1,20 +1,20 @@
 # A/B of the three automatically selected small/degenerate-shape paths against
-# the five-loop nest, via their mode overrides:
+# the five-loop nest, via the internal `path_modes` keyword of `plan_contract`:
 #
-#   outer      K = 1 outer product      C[a,b] = A[a] B[b]           _OUTER_MODE
-#   dot        M = 1 gemv               C[cde] = A[ab] B[abcde]      _DOT_MODE
-#   unpackedb  small square matmul      C[a,c] = A[a,b] B[b,c]       _UNPACKED_B_MODE
+#   outer      K = 1 outer product      C[a,b] = A[a] B[b]           outer
+#   dot        M = 1 gemv               C[cde] = A[ab] B[abcde]      dot
+#   unpackedb  small square matmul      C[a,c] = A[a,b] B[b,c]       unpacked_b
 #
 # `execute!` only (planning excluded), plus OpenBLAS `mul!` on the equivalent
 # dense matrices for scale. Modes are interleaved rep by rep; each rep times a
 # batch of calls sized to ~200 us. Prints median ns per call and [q25, q75].
 #
-#   julia --project=benchmark benchmark/bench_degenerate.jl [--reps 31]
+#   julia --project=benchmark benchmark/bench_degenerate.jl [--reps 31] [--smoke]
 
 include(joinpath(@__DIR__, "harness.jl"))
 using Statistics: quantile
 
-const REPS = argopt("reps", 31)
+const REPS = reps_arg(31)
 
 function batchsize(f!)
     f!()
@@ -39,17 +39,11 @@ end
 
 fmt((m, lo, hi)) = @sprintf("%9.0f [%6.0f,%6.0f]", m, lo, hi)
 
-function withmode(f, ref::Ref{Symbol}, mode::Symbol)
-    return () -> begin
-        old = ref[]; ref[] = mode
-        f()
-        ref[] = old
-    end
-end
-
-function row(name, ref, plan, blasf)
-    ex = () -> execute!(plan, 1.0, 0.0)
-    r = interleaved(Any[withmode(ex, ref, :never), withmode(ex, ref, :auto), blasf])
+# One plan with the path ruled out (`mode = :never`), one with it automatic.
+function row(name, mode::Symbol, args, blasf)
+    nest = plan_contract(args...; path_modes = QuasiStrided.PathModes(; mode => :never))
+    path = plan_contract(args...)
+    r = interleaved(Any[() -> execute!(nest, 1.0, 0.0), () -> execute!(path, 1.0, 0.0), blasf])
     @printf(
         "%-28s nest %s  path %s  blas %s  path/nest %.2f  path/blas %.2f\n",
         name, fmt(r[1]), fmt(r[2]), fmt(r[3]), r[2][1] / r[1][1], r[2][1] / r[3][1]
@@ -62,19 +56,19 @@ for T in (Float64, ComplexF64)
     println("## $T")
     for d in (16, 63, 128)
         a = randn(T, d); b = randn(T, d); C = zeros(T, d, d)
-        plan = plan_contract(StridedView(C), StridedView(a), (1,), StridedView(b), (2,), (1, 2))
+        args = (StridedView(C), StridedView(a), (1,), StridedView(b), (2,), (1, 2))
         am = reshape(a, d, 1); bm = reshape(b, 1, d)
-        row("outer $(d)x$(d)", QuasiStrided._OUTER_MODE, plan, () -> mul!(C, am, bm))
+        row("outer $(d)x$(d)", :outer, args, () -> mul!(C, am, bm))
     end
-    for d in (6, 8, 12, 16)
+    for d in thin((6, 8, 12, 16))
         A = randn(T, d, d); B = randn(T, d, d, d, d, d); C = zeros(T, d, d, d)
-        plan = plan_contract(StridedView(C), StridedView(A), (1, 2), StridedView(B), (1, 2, 3, 4, 5), (3, 4, 5))
+        args = (StridedView(C), StridedView(A), (1, 2), StridedView(B), (1, 2, 3, 4, 5), (3, 4, 5))
         Bm = reshape(B, d^2, d^3); av = vec(A); cv = vec(C)
-        row("dot 1x$(d^2)x$(d^3)", QuasiStrided._DOT_MODE, plan, () -> mul!(cv, transpose(Bm), av))
+        row("dot 1x$(d^2)x$(d^3)", :dot, args, () -> mul!(cv, transpose(Bm), av))
     end
-    for d in (16, 32, 64, 128)
+    for d in thin((16, 32, 64, 128))
         A = randn(T, d, d); B = randn(T, d, d); C = zeros(T, d, d)
-        plan = plan_contract(StridedView(C), StridedView(A), (1, 2), StridedView(B), (2, 3), (1, 3))
-        row("unpackedB $(d)^3", QuasiStrided._UNPACKED_B_MODE, plan, () -> mul!(C, A, B))
+        args = (StridedView(C), StridedView(A), (1, 2), StridedView(B), (2, 3), (1, 3))
+        row("unpackedB $(d)^3", :unpacked_b, args, () -> mul!(C, A, B))
     end
 end