Introduction

gemmkit is a pure-Rust workspace for GEMM, the general matrix multiply C <- alpha*A*B + beta*C. The core crate works over strided views or raw pointers. It selects the fastest instruction set available on the machine at runtime. Under a fixed input and configuration, it keeps its results reproducible from run to run. Around that core sit 3 zero-copy adapters, one each for ndarray, nalgebra, and faer. An install-time autotuner calibrates the engine for the machine it will actually run on.

This book is the narrative documentation for the whole workspace. The API reference on docs.rs remains the authority on exact signatures and item-level details. The book explains the pieces in context: the reasoning, the trade-offs, and the corners of the API that a reference page cannot cover in depth.

What is in the book

The gemmkit user guide covers the core crate. It starts with the first multiply and goes on to parts most users never need:

• matrix views and layouts
• the optional element types (f16/bf16, i8, complex)
• parallel execution
• prepacked operands
• fused epilogues
• batched GEMM
• small shapes and GEMV
• instruction-set pinning
• the tuning knobs
• no_std and WebAssembly builds
• the unchecked raw-pointer tier

The adapter guides show how to drive the engine straight from ndarray, nalgebra, and faer types with no copies. Each adapter gets a chapter with a getting-started page and an advanced page. The advanced page covers the full surface: fused operations, integer and complex GEMM, batching, and prepacking, all in the host library’s native types.

The gemmkit-tune guide explains the autotuner. It covers how to run the autotuner on a deployment machine, what the emitted profile contains, and how the sweep behind it works.

The architecture chapter walks through the inside of the engine, layer by layer. It covers how a call travels from the public API down to the microkernel, and how instruction sets and element types stay pluggable without macros. It also covers how blocking is derived from the cache hierarchy, and how the whole thing is tested. It is a more detailed, more approachable companion to the compact ARCHITECTURE.md in the repository, written to be read front to back.

How to read it

If you just want fast matrix multiplication in an application, start with Getting Started. Read the user guide as far as your use case demands.

If your matrices already live in ndarray, nalgebra, or faer, jump straight to that adapter’s chapter. When you need the underlying concepts, fall back to the user guide. The adapters forward to the same engine and share its semantics.

If you are curious how the engine works, or you plan to contribute, the architecture chapter is the intended path. It assumes you have skimmed the user guide, but not that you know BLIS. It explains the design decisions, not just the code.

Conventions and resources

Code examples target Rust edition 2024 and the workspace MSRV of 1.89. Examples that need an optional Cargo feature say so where they appear. Repository paths like gemmkit/src/driver.rs are relative to the repository root.

Related resources include the API reference, the CHANGELOG, and the crates on crates.io: gemmkit, gemmkit-ndarray, gemmkit-nalgebra, gemmkit-faer, gemmkit-tune.

本书也有简体中文版。

Getting Started

gemmkit computes C <- alpha*A*B + beta*C over strided views of ordinary Rust slices. It picks the fastest instruction set your CPU actually has, at runtime. There is no build-time ISA choice to make and no BLAS to link. You add 1 dependency, hand it 3 matrices, and call gemm.

Adding the dependency

The core crate is gemmkit. For plain f32/f64 work on a normal (std) target, this single line is all you need:

[dependencies]
gemmkit = "0.1"

This pulls in the 2 default features, std and parallel. std gives you runtime cache and CPU feature detection. It also gives you the GEMMKIT_REQUIRE_ISA and GEMMKIT_* tuning knobs, and a thread-local workspace pool. The pool makes repeated same-size calls allocation-free. parallel adds rayon multithreading, and implies std.

The optional element-type families (half, complex, int8) and the epilogue capability are off by default. A plain float build pays for none of their codegen or dependencies. To build the crate as no_std (only core + alloc), turn the defaults off with default-features = false. See no_std and WebAssembly for that setup.

A first complete example

This example computes a 2x3 matrix times a 3x2 matrix, all row-major, run single-threaded:

use gemmkit::{gemm, MatMut, MatRef, Parallelism};

fn main() {
    // 2x3 times 3x2 = 2x2, all row-major
    let a = [1.0_f32, 2.0, 3.0, 4.0, 5.0, 6.0];
    let b = [7.0_f32, 8.0, 9.0, 10.0, 11.0, 12.0];
    let mut c = [0.0_f32; 4];
    gemm(
        1.0,
        MatRef::from_row_major(&a, 2, 3),
        MatRef::from_row_major(&b, 3, 2),
        0.0,
        MatMut::from_row_major(&mut c, 2, 2),
        Parallelism::Serial,
    );
    assert_eq!(c, [58.0, 64.0, 139.0, 154.0]);
}

The arguments are exactly the terms of C <- alpha*A*B + beta*C. They are the scalar alpha, the 2 input views, the scalar beta, the output view, and a Parallelism selector.

MatRef::from_row_major(&a, 2, 3) reads a as a 2-by-3 row-major matrix. The shapes must line up: A.cols must equal B.rows, and C must be A.rows by B.cols. If they do not line up, the call panics before it touches memory.

Transposition never needs a copy. It is a stride change, not a data move. MatRef::from_col_major(&b, 3, 2) reads the same buffer as a column-major matrix. MatRef::new lets you set the row and column strides directly.

What happened under the hood

The gemm entry does a small amount of work before it does any arithmetic.

First, it validates the call. It checks that the inner dimensions agree. It checks that each view stays inside its slice. It checks that C addresses every (i, j) at a distinct offset, and that C’s storage does not overlap A’s or B’s storage. Any failed check raises a panic with a specific message, before a single unsafe operation runs. Only then does it lower the 3 views to raw pointers and strides, and hand them to the dispatch layer.

Dispatch resolves which kernel to run. The first GEMM call for a given element type runs CPU feature detection once. It records the winning entry point in a OnceLock, and returns. Every later call is a plain indirect call through that cached pointer, with no repeat detection. So the runtime ISA choice is a one-time cost, amortized across the whole process.

You can override the automatic choice with the GEMMKIT_REQUIRE_ISA environment variable. You can also use it to pin a specific backend for testing. gemmkit reads the variable once and memoizes it the same way. See Runtime ISA Dispatch for details. Life of a GEMM Call walks the full path from call to microkernel.

alpha and beta, precisely

alpha scales the product A*B. beta scales the incoming contents of C. The one detail to remember is what happens at the edges.

When beta == 0, the engine does not read C at all. It overwrites C with alpha*A*B. That rule is what makes let mut c = [0.0_f32; 4] correct above, even though you could have left the buffer uninitialized. In concrete terms, a beta == 0 output slice may hold garbage. Through the unchecked tier it may even be genuinely uninitialized memory, and the result is still well defined.

When beta == 1, the engine leaves the existing C untouched and accumulates the product onto it. Any other beta value first multiplies C through.

There is also a degenerate fast path. If k == 0 (an empty contraction) or alpha == 0 (the product vanishes), the call reduces to C <- beta*C. It never touches A or B at all, and just scales the output in place. Combined with the beta == 0 rule, alpha == 0, beta == 0 zeroes C, and k == 0, beta == 1 is a no-op. Narrow types scale in f32 and round back on the store. So the degenerate path rounds exactly as the full kernel would.

The Cargo features

The element-type and capability features compose. half + epilogue gives fused f16 GEMM. int8 + epilogue gives the requantizing entries, and so on. See Element Types and Fused Epilogues for each combination.

Version requirements

gemmkit targets Rust 1.89 on edition 2024. It is licensed MIT OR Apache-2.0. The API reference is on docs.rs/gemmkit. This book is the long-form companion to that reference.

Where to next

• Matrix Views and Layouts: how MatRef/MatMut, strides, transposition, and submatrices work, and exactly what the safe API validates.
• Element Types: f16/bf16, i8, and complex, with their accuracy characteristics.
• Parallelism in Practice: what Rayon(0) auto actually does and when Serial is the right call.
• Runtime ISA Dispatch and The Unchecked Tier: overriding the backend, and the raw-pointer engine.
• Working with an existing array library? See the adapters: ndarray, nalgebra, faer.

Matrix Views and Layouts

Every gemmkit call takes its operands as views: a slice, a shape, and 2 strides. MatRef<'a, T> is the immutable input view. MatMut<'a, T> is the mutable output view. Neither view owns its data. Both borrow a slice you already have.

The library’s whole layout vocabulary lives in those 2 stride numbers: row-major, column-major, transposed, submatrix, and broadcast. So the same buffer can be read a dozen ways without ever being copied.

The 2 strides

Element (i, j) of a view lives at slice offset i*rs + j*cs, where rs is the row stride and cs is the column stride. Strides are counted in elements, not bytes. An rs of 4 means the next row sits 4 elements further along the slice. That single offset formula is the entire model. Everything else is a choice of rs and cs.

3 constructors cover the common cases. Each exists on both MatRef and MatMut:

#![allow(unused)]
fn main() {
use gemmkit::MatRef;

let data = [0.0_f32; 12];
let row_major = MatRef::from_row_major(&data, 3, 4); // rs = cols = 4, cs = 1
let col_major = MatRef::from_col_major(&data, 3, 4); // rs = 1, cs = rows = 3
let general   = MatRef::new(&data, 3, 4, 4, 1);      // explicit rs, cs (here == row-major)
}

from_row_major(data, rows, cols) sets rs = cols, cs = 1. Rows are contiguous: the classic C order. from_col_major(data, rows, cols) sets rs = 1, cs = rows. Columns are contiguous: Fortran order.

new(data, rows, cols, rs, cs) takes the strides verbatim. Reach for it when neither canonical layout matches, for example a submatrix, or a view whose leading dimension differs from its logical width. MatRef and MatMut also expose .rows() and .cols().

Transposition is a stride swap

Because (i, j) maps through i*rs + j*cs, swapping the roles of the 2 strides (and the 2 dimensions) transposes the view in place. Say a holds an m x k matrix in row-major order (rs = k, cs = 1). Its transpose is the k x m matrix whose (i, j) is the original (j, i), at offset j*k + i. That offset is exactly rs = 1, cs = k over the same slice:

#![allow(unused)]
fn main() {
use gemmkit::MatRef;

// `a` is m x k row-major
let (m, k) = (2, 3);
let a = [1.0_f32, 2.0, 3.0, 4.0, 5.0, 6.0];

let a_rowmajor = MatRef::from_row_major(&a, m, k); // m x k
let a_transposed = MatRef::from_col_major(&a, k, m); // k x m, same bytes, no copy
}

So from_col_major over a row-major buffer is the transpose, and the reverse also holds. new with rs/cs swapped does the same for any layout. A transposed operand therefore costs nothing at the API level: the kernel walks the strides you give it. This is how you feed A^T * B or A * B^T without materializing a transpose.

Submatrices and strided views

A submatrix is a view whose leading dimension (the distance between successive rows or columns) is larger than its logical extent. Build one by slicing the buffer so the block’s top-left element sits at the start of the slice, then handing over the parent’s strides. Here is the top-left 2 x 2 block of a 4 x 4 row-major matrix, starting at row 1, column 1:

#![allow(unused)]
fn main() {
use gemmkit::MatRef;

let parent = [0.0_f32; 16]; // 4x4 row-major, leading dimension 4
let block = MatRef::new(&parent[1 * 4 + 1..], 2, 2, 4, 1); // rs stays 4, cs stays 1
}

The row stride is still 4, the parent’s width. So consecutive rows of the block skip over the columns you excluded. The slice begins at offset 5, the block’s (0, 0). The safe API verifies the tail slice is long enough to reach the block’s far corner.

The same mechanism expresses a broadcast input. A stride of 0 makes a dimension repeat 1 element. A 1 x n row broadcast down m rows is MatRef::new(row, m, n, 0, 1). Every logical row then reads the same storage. gemmkit allows broadcasts for the read-only inputs A and B, but never for the output. The next section explains why.

What the safe API accepts, and what it rejects

The safe entries (gemm, gemm_i8, gemm_cplx, and the fused variants) accept non-negative strides only, including 0 for a broadcast input. A negative stride is outside what a &[T] view can describe safely. So is a base pointer that sits in the middle of a buffer, rather than at element (0, 0). Those cases live in The Unchecked Tier, the raw-pointer engine the adapters use to express arbitrary layouts.

Before any arithmetic, the safe entries run one validation prologue over the (A, B, C) trio. Every failure is a panic, raised ahead of the first unsafe operation:

• Shape agreement. A.cols == B.rows, A.rows == C.rows, B.cols == C.cols. A mismatch panics with the offending pair, for example gemmkit: A.cols (3) != B.rows (4).
• In-bounds views. For each view, the engine computes the highest slice offset it will touch, and checks it against the slice length. Too small a slice panics with gemmkit: A view of 3x4 (strides 4,1) needs 12 elements but slice has 8. A view whose strides are negative, or so large the addressing overflows usize, panics with ... has negative strides or is too large to address; use gemm_unchecked.
• C addresses each element uniquely. gemmkit writes the output, so 2 distinct (i, j) must never land on the same offset. A self-aliasing C, such as a zero row or column stride, or strides that collide, would become a data race in parallel mode. That is reachable from entirely safe code, so it panics: gemmkit: C view aliases itself (...); C must address each (i,j) uniquely. This is why broadcast strides are fine for A/B (read-only) but forbidden for C.
• C does not overlap A or B. The output’s byte range must be disjoint from each input’s. gemmkit compares byte ranges, not element counts, so the heterogeneous integer API (i8 inputs, i32 output) stays exact. Overlap panics with gemmkit: C aliases A or B. In fully safe Rust, the borrow checker already forbids an overlapping &mut/& pair. This check is a defensive backstop that also covers the raw-lowered paths.

These messages are stable. The correctness suite asserts on their wording, so you can rely on them in tests.

Zero-sized dimensions

A view with a zero dimension is legal and validates cleanly. gemmkit accepts a 0 x k, m x 0, or m x n x (k = 0) shape. Any slice, even an empty one, satisfies the in-bounds check, because such a view addresses nothing.

If m == 0 or n == 0, the call is a no-op: there is no output to write. If only k == 0, the contraction is empty, and the call reduces to C <- beta*C. That is the same scale-only path alpha == 0 takes. See Getting Started for that degenerate rule.

Where to next

• Element Types: the same views over f16/bf16, i8, and complex data.
• The Unchecked Tier: negative strides, interior base pointers, and the raw-pointer engine.
• The adapters (ndarray, nalgebra, faer) build these views for you from each library’s native matrix types.

Element Types

gemmkit multiplies more than f32. The same engine, driver, and blocking model serve 4 element-type families. Each family is a Cargo feature. Each also has a SIMD implementation on every backend, over the portable scalar fallback.

What changes between the families is the input type, the accumulator type, and the output type. With that comes a change in the accuracy you should expect. This page maps what is available and how precise each family is.

The built-in real floats

f32 and f64 need no feature flag. They go through the generic gemm (and gemm_with, and the unchecked entries), and accumulate in their own type. They are the baseline every other family is measured against:

#![allow(unused)]
fn main() {
use gemmkit::{gemm, MatMut, MatRef, Parallelism};

let a = [1.0_f64, 2.0, 3.0, 4.0];
let b = [5.0_f64, 6.0, 7.0, 8.0];
let mut c = [0.0_f64; 4];
gemm(2.0, MatRef::from_row_major(&a, 2, 2), MatRef::from_row_major(&b, 2, 2),
     0.0, MatMut::from_row_major(&mut c, 2, 2), Parallelism::Serial);
}

Accuracy follows the textbook GEMM story. Relative error grows roughly with the contraction depth k and the machine epsilon of the type. The correctness suite holds results to a relative Frobenius gate of 8*k*eps, checked against an independent f64 reference. So f64 is near exact for any realistic k, and f32 carries its usual ~1e-7 per-element relative precision.

Narrow floats: the half feature

With half on, f16 and bf16 become element types. gemmkit re-exports them as gemmkit::f16 and gemmkit::bf16, so you need not depend on half directly. They share the generic gemm surface: MatRef<'_, f16> in, MatMut<'_, f16> out, because they implement the same scalar trait as the real floats.

The defining property is mixed precision. The engine widens inputs to f32 on load. The entire contraction accumulates in f32. The engine rounds the result back to the narrow type exactly once, at the store. There is no repeated narrow rounding inside the k loop. That single rounding point is what keeps the accuracy usable.

#![allow(unused)]
fn main() {
use gemmkit::{f16, gemm, MatMut, MatRef, Parallelism};

let a: Vec<f16> = (0..6).map(|i| f16::from_f32(i as f32)).collect();
let b: Vec<f16> = (0..6).map(|i| f16::from_f32(i as f32)).collect();
let mut c = vec![f16::ZERO; 4];
gemm(f16::ONE, MatRef::from_row_major(&a, 2, 3), MatRef::from_row_major(&b, 3, 2),
     f16::ZERO, MatMut::from_row_major(&mut c, 2, 2), Parallelism::Serial);
}

Because the accumulation is in f32, the dominant error is that single final round, not the sum. f16 carries about 9.8e-4 (2^-10) relative precision, and bf16 about 7.8e-3 (2^-7), both essentially independent of k. A narrow-precision GEMM is therefore close to computing in f32 and rounding once. That is far more accurate than accumulating in 16 bits would be.

One consequence of rounding once is that at a large k, a single depth panel would stream an intermediate result too large for L2 cache. The engine handles this itself. Past an auto-derived byte gate, it switches to an f32-output internal twin. That twin re-blocks the contraction to stay cache-resident, and narrows the result at the end. The twin matches the single panel byte for byte for the common case of beta in {0, 1}, and stays within tolerance otherwise. This switch is automatic and needs no configuration. The mechanism is detailed in Dot Kernels and the Deep-K Twin. On AVX-512 BF16 hardware, bf16 also uses the vdpbf16ps dot kernel. See Runtime ISA Dispatch for that mechanism.

Integer: the int8 feature

int8 adds gemm_i8, a separate entry point. Its input and output types differ, i8 in and i32 out, and the homogeneous gemm<T> surface cannot express that. alpha, beta, and C are all i32:

#![allow(unused)]
fn main() {
use gemmkit::{gemm_i8, MatMut, MatRef, Parallelism};

let a = [1_i8, 2, 3, 4, 5, 6];
let b = [7_i8, 8, 9, 10, 11, 12];
let mut c = [0_i32; 4];
gemm_i8(1, MatRef::from_row_major(&a, 2, 3), MatRef::from_row_major(&b, 3, 2),
        0, MatMut::from_row_major(&mut c, 2, 2), Parallelism::Serial);
}

Integer GEMM is exact. It is i32 ring arithmetic that wraps on overflow, the conventional integer-GEMM semantics. There is no tolerance to speak of, because there is no rounding. The result is bit-for-bit identical across every ISA (scalar, FMA, AVX-512F, and the AVX-512 VNNI vpdpbusd dot kernel). It is also identical between a serial and a parallel run, because integer addition over a ring does not depend on order.

If you feed values whose products can exceed i32, the wraparound is defined and reproducible. It is not undefined behavior. The int8 feature pulls in no extra dependency. Adding epilogue on top unlocks the requantizing entries, which give i8 or u8 output in one pass. See Fused Epilogues for those entries.

Complex: the complex feature

complex adds gemm_cplx over num-complex values. gemmkit re-exports them as gemmkit::c32 (Complex<f32>) and gemmkit::c64 (Complex<f64>). Its signature carries a conjugation flag for each operand:

#![allow(unused)]
fn main() {
use gemmkit::{c32, gemm_cplx, Complex, MatMut, MatRef, Parallelism};

let a = [Complex::new(1.0_f32, 1.0), Complex::new(2.0, 0.0)];
let b = [Complex::new(0.0_f32, 1.0), Complex::new(1.0, 0.0)];
let mut c = [c32::default(); 1];
gemm_cplx(
    Complex::new(1.0, 0.0),
    MatRef::from_row_major(&a, 1, 2), false, // conj_a
    MatRef::from_row_major(&b, 2, 1), false, // conj_b
    Complex::new(0.0, 0.0),
    MatMut::from_row_major(&mut c, 1, 1),
    Parallelism::Serial,
);
}

The computation is C <- alpha*op(A)*op(B) + beta*C. op(A) is conj(A) when conj_a is set, and likewise op(B) is conj(B) when conj_b is set. Passing false, false gives the plain product A*B. The flags conjugate the operands only.

Complex accumulates in its own type. It is held to a relative Frobenius gate of 16*k*eps, with eps the real component’s epsilon. So a c32 GEMM is about as accurate as an f32 one, and a c64 GEMM about as accurate as an f64 one. complex pulls in num-complex.

Internally, complex does not ride the float kernel. It uses a dedicated split (structure-of-arrays) kernel instead, which is why it gets a separate entry point. That design is covered in The Complex Split Kernel.

Choosing a type

If you need speed and can tolerate ~1e-3 precision, choose bf16 or f16. Both use half the bytes per element that f32 does, and both still accumulate in f32. If you need exactness, int8 gives it. If you need range and precision, stay on f32/f64.

The reproducibility contract is the same for all of them. See Parallelism in Practice for what “reproducible” does and does not promise.

Where to next

• Fused Epilogues: bias, activation, requantization, and per-element maps fused into the store.
• Runtime ISA Dispatch: the VNNI and BF16 dot kernels these types can reach.
• Matrix Views and Layouts: constructing the views every family shares.

Parallelism in Practice

Every GEMM entry takes a Parallelism argument as its last parameter. It is a small enum with 3 practical modes. Using it well comes down to 2 things: understand what the auto mode decides for you, and know when to take over yourself.

The 3 modes

#![allow(unused)]
fn main() {
pub enum Parallelism {
    Serial,       // single-threaded
    Rayon(usize), // rayon with at most n threads; Rayon(0) auto-detects
}
}

Serial runs the whole call on the calling thread. Rayon(n) asks for at most n workers. Rayon(0) is auto, and it is also the Default, so Parallelism::default() gives you auto. The n in Rayon(n) is a ceiling on partitions, not a promise to use them all. A problem with less work than n chunks, or fewer cores than n, gets fewer.

What auto actually does

Auto does not mean “use all cores.” It makes 2 decisions based on the problem size.

First, a workload gate applies. Below a total-work threshold on m*n*k (the GEMMKIT_PARALLEL_THRESHOLD knob, default 48*48*256), the call stays serial no matter what. On a matrix that small, fork/join overhead would swamp any gain. This gate runs before everything else, so it applies even to an explicit Rayon(n). Below the gate, Rayon(8) still runs on one thread.

Above the gate, auto scales the worker count with the total work instead of jumping straight to the full core count. It targets m*n*k divided by GEMMKIT_PAR_MNK_PER_WORKER (default 2_000_000, one worker per that much work). It then caps the result by the machine’s core count and by the number of available job chunks, floored at one. The count is work-based, not dimension-based, because the best worker count tracks total flops, not linear size. No single stride on a linear dimension can fit that whole range. A small product uses a handful of workers, and a large one uses many. Setting GEMMKIT_PAR_MNK_PER_WORKER to 0 (which behaves as 1) forces full width for anything above the serial gate.

Explicit counts

Rayon(n) with n > 0 bypasses the ramp heuristic and asks for exactly n partitions. For safety, this is still capped by the machine’s core count (available_parallelism) and by the number of job chunks the problem actually splits into. So Rayon(1000) on a 16-core box computing a small product does not oversubscribe: it collapses to what the machine and the work can absorb. This exactness is why the test suite and the scaling diagnostics use explicit counts. Rayon(4) gives you four-way partitioning, when there is that much work and that many cores, not a heuristic guess. Use an explicit count once you have measured your own workload and know the sweet spot. Also use it when you want reproducible partitioning across runs for benchmarking.

How gemmkit uses the rayon pool

gemmkit does not require you to hand it a rayon pool. If you wrap a call in your own pool’s install, the GEMM’s workers run on that pool instead of anywhere else.

#![allow(unused)]
fn main() {
let pool = rayon::ThreadPoolBuilder::new().num_threads(4).build().unwrap();
pool.install(|| {
    gemm(1.0, a, b, 0.0, c, Parallelism::Rayon(0)); // runs on `pool`
});
}

The worker count gemmkit chooses is still bounded by available_parallelism, the whole machine. Rayon’s work-stealing scheduler distributes the partitions over whatever threads the current pool has. A smaller custom pool simply runs the same partitions on fewer threads. Work distribution inside a call is demand-driven: workers pull contiguous chunks from a shared lock-free cursor. On a heterogeneous part (a mix of P-cores and E-cores), a faster core absorbs more chunks instead of everyone waiting on the slowest one.

If a call does not run inside a pool you installed yourself, gemmkit reaches for one of its own pools instead. This is on by default on native targets, x86_64 and aarch64. gemmkit keeps up to GEMMKIT_POOL_CLASSES (default 2 on x86_64, 1 on aarch64) private, persistent pools. Each is sized to an exact halving tier of the machine width. On a 32-thread part that means tiers of 16 and 8 threads. On a 14-core M4 Max it means a single 7-wide tier. Each pool builds lazily on first use and is never rebuilt.

Auto snaps its worker count exactly onto one of these tiers instead of forking the full-width global pool. A fork’s overhead tracks the pool’s idle slack: the threads it owns beyond the ones actually engaged. A small GEMM drowns in that slack on a full-width pool, so matching the pool to the work avoids the drag.

None of this changes what you already know above. An install’d call is still fully respected and never redirected to a tier pool. An explicit Rayon(n) still gets exactly n workers, just routed into whichever tier pool is the smallest fit. What does change is idle memory. By default, an x86_64 process now parks about 24 extra threads, the 16- and 8-wide tier pools, alongside the global pool. An aarch64 M4 Max process parks 7 (its single half-width tier). All of them stay asleep until a small GEMM needs them. Set GEMMKIT_POOL_CLASSES=0 to disable tier pools entirely. Every call then falls back to the ambient pool.

The threaded-wasm story is different: there gemmkit always sizes a dedicated pool of its own. See no_std and WebAssembly for that case.

The reproducibility promise, precisely

For a fixed input, environment, and configuration, the output is identical regardless of the worker count. That is the contract.

It holds for 2 reasons. First, kc, nc, and the fixed depth-panel order are the only things that shape each output element’s summation. gemmkit computes all 3 independently of how many threads will run them. Second, a single worker reduces each output element start to finish, over the full contraction depth. No split reduction exists whose order could depend on the schedule.

The flat job list itself is not strictly identical across worker counts: a wide worker count can shrink mc to keep the list deep enough. But mc always stays a multiple of mr, so the set of microtiles, and their numerics, stay unchanged. The packed bytes do not depend on who packs them, either. Which worker computes a given tile varies from run to run. The numerical result does not.

What is not promised is bitwise identity between Serial and Rayon(n). It happens to hold today on the driver paths, since serial and parallel run the same kernel. But build on reproducibility under a fixed config, not on serial-versus-parallel bit equality. You will not get cross-machine or cross-config bit equality here. Floating-point GEMM is order-sensitive, and the config (ISA, blocking, thread cap) is part of the fixed input. Integer gemm_i8 is the exception. It is bit-identical across ISAs and worker counts, because i32 addition is order-independent.

When Serial is the right call

Reach for Serial in 3 situations.

1. Small problems. Below the workload gate, auto is serial anyway. Passing Serial explicitly also skips the available_parallelism probe and the fork machinery entirely, which is cheaper in a tight loop of tiny GEMMs.
2. When you own the outer parallelism. Suppose you already run many independent GEMMs across a rayon pool, or you parallelize a batch loop yourself. Do not let each inner call also fan out. That oversubscribes the machine, and it usually hurts performance instead of helping. Run the inner calls Serial and keep the parallelism at the outer level. For a batch of products, prefer the built-in Batched GEMM entries instead, since they schedule the whole batch as one unit.
3. Determinism-sensitive debugging. Use the single-threaded path to rule out scheduling as a variable.

Bandwidth-bound shapes get their own policy

A matrix-vector product (m == 1 or n == 1), and other memory-bound shapes, are not compute-bound. The work-based worker count above is the wrong model for them, so these routes use a separate rule.

Below a cache-derived byte floor, the matrix fits one core’s private cache. That core already saturates the cache alone, so splitting the work only adds contention, and the route stays serial. Above the floor, the route steps straight to a width chosen for the bytes touched. This width tops out at half the logical core count, because a gemv saturates its bandwidth well before the last core joins in.

That width climbs in steps, not smoothly. Each step is one of the exact-fit thread pools described above, so a bandwidth-bound call also gets a pool sized exactly to it. A few workers is the worst point on a bandwidth scaling curve, so the policy jumps over that point instead of ramping through it.

This whole policy is automatic. The byte floor, the step spacing, and a flat override are all tunable, through GEMMKIT_GEMV_PARALLEL_BYTES, GEMMKIT_GEMV_TIER_STEP, and GEMMKIT_GEMV_THREAD_CAP. See Small Shapes and GEMV for the full treatment.

Where to next

• Small Shapes and GEMV: the bandwidth-bound policy in detail.
• Batched GEMM: scheduling many products as one batch instead of nesting parallelism.
• Tuning Knobs: the GEMMKIT_* thresholds behind these decisions.
• Parallel Execution: the job cursor and worker-count resolution internals.

Prepacked Operands

Before the microkernel can touch A and B, the engine copies each into a cache-friendly micropanel layout. This layout holds contiguous tiles, and the microkernel walks each tile with unit strides. For a one-shot product, that copy is pure setup. The engine pays for it once and never uses it again.

Many workloads multiply the same matrix over and over. A linear layer applies one fixed weight matrix to a stream of activation batches. A solver applies one fixed operator to many right-hand sides. Repacking a fixed operand on every call throws away work the engine already did. The prepacked-operand API lets you pay for the pack once, then reuse the result across every product that shares that operand.

Packing the right-hand side

The common case fixes B (the weights) and streams a series of differently sized A matrices (the activations). Call prepack_rhs once to turn a k x n B into a PackedRhs handle. Then feed that handle to gemm_packed_b for each product:

#![allow(unused)]
fn main() {
use gemmkit::{prepack_rhs, gemm_packed_b, MatRef, MatMut, Parallelism};

// fixed weights: a k x n matrix reused across many activation batches
let (k, n) = (512, 256);
let weights = vec![0.0f32; k * n];
let packed = prepack_rhs(MatRef::from_col_major(&weights, k, n));

// per activation batch: an m x k input, sharing the packed weights
let mut c = vec![0.0f32; m * n];
gemm_packed_b(
    1.0,
    MatRef::from_row_major(&input, m, k),
    &packed,
    0.0,
    MatMut::from_col_major(&mut c, m, n),
    Parallelism::Rayon(0),
);
}

prepack_rhs accepts any layout of B and reads it through its strides. A row-major, column-major, or transposed view all pack the same way. The pack runs once, inside prepack_rhs, on a single thread. Every later call to gemm_packed_b skips it.

The buffer records its own blocking geometry: nr, kc, and nc. Every consuming call reads that same geometry back, so a panel is always read against its own tiling. The buffer is read-only for the whole GEMM, so gemmkit never writes to it after the pack. A single PackedRhs is therefore safe to share across threads and across concurrent calls, with no extra synchronization. PackedRhs::rows() reports the original k. PackedRhs::cols() reports the original n.

The handle works for any product whose (k, n) match the packed B, as long as the output C is column-major-ish (|csc| >= |rsc|). That constraint is the one surprise in this API. A row-major C would force the engine to swap A and B internally, to keep its stores contiguous. A prepacked B cannot move into the A role. So gemm_packed_b panics on a row-major C and points you back to plain gemm for that layout. Only C is pinned this way. A’s layout is unconstrained.

Under a fixed configuration, gemm_packed_b reproduces plain gemm and is deterministic across worker counts. There is one narrow caveat. It applies to a small product, where both m and n are at or below the small_mn_dim knob (16 by default, 32 on aarch64). It also applies to a gemv-shaped product, where m == 1 or n == 1. In both cases the 2 calls may differ in the last ULP. The reason is routing, not error. Plain gemm reroutes those shapes to a special path, while the prepacked entry always drives the general packed kernel. Both answers are correct. They sum in a slightly different order, only on the shapes where the special paths would otherwise take over.

The left-hand-side mirror

The symmetric case fixes A and streams a series of varying B matrices. prepack_lhs produces a PackedLhs handle, and gemm_packed_a consumes it. This mirrors the RHS pair exactly, with the axes relabeled. PackedLhs::rows() is the original m. PackedLhs::cols() is the shared k.

Internally, the LHS pack is not a separate code path. By the engine’s A/B symmetry, a prepacked A is exactly the prepacked B of the transposed product C^T = B^T A^T. So prepack_lhs lays down the identical micropanel buffer, and only records the dimensions in LHS terms.

This has one visible consequence: the C-layout constraint flips. gemm_packed_a requires a row-major-ish C (|csc| <= |rsc|), the exact opposite of the RHS entry. A column-major C would keep A in the genuine LHS role, and a transposed-RHS buffer cannot fill that role. Pick the packed-A entry when your C is row-major. Pick the packed-B entry when your C is column-major. Together the 2 entries cover both orientations.

Fused variants

Each packed entry has a fused twin, under the epilogue feature. gemm_packed_b_fused and gemm_packed_a_fused add a per-row or per-col bias, plus an optional activation, in the same store the packed kernel already runs. See Fused Epilogues for the bias and activation types.

The same PackedRhs or PackedLhs handle serves both the plain entry and the fused entry. The epilogue applies only at the store, so it never touches the pack. Build the buffer once, then choose per call whether to fuse.

2 details are specific to the packed path. First, unlike plain gemm_fused, the packed fused entries never reroute to the gemv, small-m,n, or small-k kernels. They always drive the general packed kernel, the same divergence the plain packed entries document. Second, gemmkit always gives the per-row or per-col bias in the natural user frame. gemm_packed_a_fused handles the internal transpose for you, so a PerRow bias has length A.rows, no matter which entry you call.

Prepacking i8 weights

Under the int8 feature, the same pattern extends to quantized inference. prepack_rhs_i8 packs a fixed i8 weight matrix into a PackedRhs<i8>, and gemm_i8_packed_b consumes it. It takes i8 inputs and produces an i32 output.

Prepacking is a bigger win here than for floats, for a structural reason. The AVX-512 VNNI kernel (vpdpbusd) reads its RHS from a k-quad-interleaved layout. The engine cannot produce that layout in place, so this kernel’s RHS pack is mandatory on every call. At small m, that per-call O(k*n) pack easily dominates the O(m*k*n) compute. Prepacking removes it from the hot loop entirely.

The packed buffer also pins the kernel choice. It is laid out for whichever integer kernel the process’s dispatch selected, either the VNNI interleave or the widen kernel’s plain panels. gemm_i8_packed_b always runs that same family, so the buffer is never misread.

Integer accumulation is exact and does not depend on the ISA. So the packed and plain paths agree bit-for-bit for every valid shape, with no small-shape caveat at all.

#![allow(unused)]
fn main() {
use gemmkit::{prepack_rhs_i8, gemm_i8_packed_b, MatRef, MatMut, Parallelism};

let packed = prepack_rhs_i8(MatRef::from_col_major(&weights_i8, k, n));
let mut c = vec![0i32; m * n];
gemm_i8_packed_b(
    1,
    MatRef::from_row_major(&input_i8, m, k),
    &packed,
    0,
    MatMut::from_col_major(&mut c, m, n),
    Parallelism::Rayon(0),
);
}

When prepacking pays

Prepacking trades one upfront O(k*n) copy for a saved repack on every later product against that operand. It pays off exactly when you reuse the operand enough to amortize that copy. A single product, or an operand that changes on every call, gains nothing. The pack is then pure overhead, and plain gemm is the right tool.

Be aware that the float path does not always pack B in the first place. For small m, plain gemm reads B in place, unpacked, a choice governed by the rhs_pack_threshold knob. So prepacking a lightly reused float B can even lose.

The clearest wins are the fixed-weight inference loops this API is named for. Above all, that means the i8 VNNI path, whose RHS pack is otherwise unavoidable on every single call. When in doubt, measure the loop both ways. The crossover point depends on your reuse count and your machine, not on a fixed rule.

The raw-pointer forms (prepack_rhs_unchecked, gemm_packed_b_unchecked, and their _with, LHS, and i8 counterparts) exist for adapters and FFI that validate their own inputs. See The Unchecked Tier.

Fused Epilogues

A GEMM rarely stands alone. Its output usually feeds straight into a bias add, an activation, or a quantization step. Done naively, that means a second full pass over C. The GEMM writes m*n values, then a separate loop reads them all back, transforms them, and writes them again.

A fused epilogue folds that transform into the GEMM’s own store. It transforms each output element in-register, at the moment the store writes it, so the extra pass over memory disappears. Everything on this page lives behind the epilogue Cargo feature.

Bias and activation

gemm_fused is the vectorized workhorse. It computes C <- act(alpha*A*B + beta*C + bias) in one pass. The bias is a Bias enum: either Bias::PerRow(&[T]) (one value per output row, length m) or Bias::PerCol(&[T]) (one per output column, length n). gemmkit adds that value to every element of the matching row or column, after the product. The activation is an Activation: Relu (max(v, 0)) or LeakyRelu(slope). Both arguments are Option. Passing None for both delegates straight to plain gemm.

#![allow(unused)]
fn main() {
use gemmkit::{gemm_fused, Bias, Activation, MatRef, MatMut, Parallelism};

let bias = vec![0.0f32; m]; // one value per output row
gemm_fused(
    1.0,
    MatRef::from_row_major(&a, m, k),
    MatRef::from_col_major(&b, k, n),
    0.0,
    MatMut::from_col_major(&mut c, m, n),
    Some(Bias::PerRow(&bias)),
    Some(Activation::Relu),
    Parallelism::Rayon(0),
);
}

The bias, the LeakyRelu slope, and the activation all apply in-register on the vector fast path. So the fusion costs almost nothing over the raw GEMM.

An arbitrary per-element map

When the transform is not a bias or a standard activation, use gemm_map. It takes a closure f(value, row, col) -> value. gemmkit applies that closure to each output element, exactly once, at its final value, fused into the store. It is the general extension point for an epilogue with no dedicated fast path in gemmkit: GELU, sigmoid, a clamp, or a position-dependent transform.

#![allow(unused)]
fn main() {
use gemmkit::{gemm_map, MatRef, MatMut, Parallelism};

let f = |v: f32, _r: usize, _c: usize| v.tanh();
gemm_map(
    1.0,
    MatRef::from_row_major(&a, m, k),
    MatRef::from_col_major(&b, k, n),
    0.0,
    MatMut::from_col_major(&mut c, m, n),
    &f,
    Parallelism::Rayon(0),
);
}

The (row, col) handed to the closure are in the user frame of C. The closure may capture its environment by reference. The bound is + Sync, so gemmkit can share that reference safely across the parallel workers, for example to borrow a lookup table. gemm_map works for f32/f64 only. It trades one indirect call per output element, cheap next to the O(k) work per element, for total generality. For a plain bias or activation, prefer gemm_fused, which vectorizes the transform.

Integer requantization

Quantized inference wants the opposite of a widening GEMM. It takes i8 inputs and accumulates into i32. It produces an i8 (or u8) output again, applying a scale and a zero-point on the way down. gemm_i8_requant and gemm_i8_requant_u8 do the whole thing in one pass. That deletes the full m*n i32 materialization a separate gemm_i8 call, followed by a requantize step, would need. Both entries take a Requantize struct:

#![allow(unused)]
fn main() {
use gemmkit::{gemm_i8_requant_u8, Requantize, RequantScale, MatRef, MatMut, Parallelism};

let req = Requantize {
    scale: RequantScale::PerRow(&per_channel_scales), // length m, per-channel
    zero_point: 128,
    bias: Some(&i32_bias),                             // optional per-row i32 bias, length m
};
gemm_i8_requant_u8(
    MatRef::from_row_major(&activations, m, k),
    MatRef::from_col_major(&weights, k, n),
    req,
    MatMut::from_col_major(&mut out_u8, m, n),
    Parallelism::Rayon(0),
);
}

The output is C[i,j] = clamp(zero_point + round_ne(scale * (sum_k A*B + bias[i])), LO, HI), using round-half-to-even. scale is either a single RequantScale::PerTensor(f32) or a per-row RequantScale::PerRow(&[f32]), the per-channel convention. The entry sets the clamp band: [-128, 127] for gemm_i8_requant, [0, 255] for the u8 twin. There is no alpha, since it folds into scale. There is no beta, since accumulating into an already-quantized C is not well defined. The requantize map is bit-exact across every ISA (scalar, FMA, AVX-512F, VNNI) and across the vector and scalar store paths. So the answer never depends on which kernel ran.

Complex bias

Under the complex feature, gemm_cplx_fused adds a per-row or per-col bias to a complex product: C <- alpha*op(A)*op(B) + beta*C + bias. It takes the same optional operand conjugation as gemm_cplx. It is bias-only by design. An ordering-based activation like ReLU has no definition on complex numbers. The conj_a and conj_b flags conjugate the operands only. gemmkit adds the bias verbatim and never conjugates it.

What you can rely on

Every fused entry routes each shape through the same kernel plain gemm would pick: the general driver, or one of the special paths. It fuses the epilogue into that kernel’s store without changing its accumulation order. So a fused call is not a different algorithm. It runs the same GEMM and applies the map at store time. The concrete guarantees:

• For f32/f64, a fused result is bitwise identical to plain gemm followed by the same scalar map. This holds for every shape, every layout, and every worker count. gemm_map gives the same guarantee for a per-element f. The complex bias entry gives the same guarantee against a complex gemm_cplx call followed by the bias add.
• For the narrow floats f16/bf16 (feature half), there is one documented exception. gemmkit widens the bias and the slope exactly to f32, applies the epilogue in f32, and narrows to the output once, on store, with round-to-nearest-even. That is more precise than gemm followed by a separate map, which would round to the narrow type, widen, and round again. So for narrow types, the fused result is deliberately not bitwise-equal to that two-step form. Reproducibility and determinism are unchanged.
• Serial and parallel runs agree bit-for-bit today. The identity-fused case (None/None, or an absent bias) const-folds back to exactly plain gemm. The reproducibility contract covers only a fixed configuration, and the worker count is part of that configuration.

The payoff is the pass over C you no longer make. On a memory-bound epilogue, that second pass can cost as much as the store itself. So fusing a bias or activation into the GEMM is close to free, in a case where the two-step form is not.

The fused epilogues also compose with the other API tiers. gemm_batched_fused applies one shared bias and activation to every element of a batched GEMM. gemm_packed_b_fused and gemm_packed_a_fused fuse over a prepacked operand. Every checked entry has raw-pointer _unchecked twins for adapters and FFI. Those twins carry the bias as a (ptr, BiasDim) pair instead of the Bias enum. See The Unchecked Tier.

Batched GEMM

A single large GEMM saturates a modern CPU on its own. A crowd of tiny ones does not. Attention heads, grouped convolutions, per-sample linear layers, and block-diagonal solves all produce many small independent products. Running them as a plain loop of gemm calls leaves most of the machine idle. Each call is too small to parallelize usefully. Yet the loop still pays a fork/join, or serializes, once per element. The batched entries take the whole set in one call and schedule it as a unit. They assign whole GEMMs to workers, so a batch of small matrices actually fills the cores.

Batched GEMM is an orchestration layer, not a new kernel. Each element re-dispatches through the full single-GEMM engine. So a batch composes automatically with the driver, the gemv path, and the small-shape paths. A batch of 1 x 1 x k products runs the horizontal dot inside each element. A batch of ordinary shapes runs the register-blocked driver. Every element is an independent GEMM, and the whole batch is reproducible across worker counts.

The strided form

The elements might sit at a regular stride, one matrix after another in a flat buffer. In that case, gemm_batched takes the single-element shape and strides once, plus a batch stride for each of A, B, and C. Element b sits at A + b*a_batch_stride, B + b*b_batch_stride, and C + b*c_batch_stride. All elements share the same shape:

#![allow(unused)]
fn main() {
use gemmkit::{gemm_batched, MatRef, MatMut, Parallelism};

// `batch` independent m x k times k x n products, packed contiguously
gemm_batched(
    batch,
    1.0,
    MatRef::new(&a, m, k, 1, m as isize), (m * k) as isize, // A element + batch stride
    MatRef::new(&b, k, n, 1, k as isize), (k * n) as isize, // B element + batch stride
    0.0,
    MatMut::new(&mut c, m, n, 1, m as isize), (m * n) as isize, // C element + batch stride
    Parallelism::Rayon(0),
);
}

A batch stride of 0 broadcasts one operand across the whole batch. This is valid for the read-only A or B, for example one shared weight matrix against a batch of inputs, but never for C. Workers write C’s elements concurrently, so those elements must stay disjoint. The result reproduces a loop of gemm calls exactly.

Under the epilogue feature, gemm_batched_fused applies one shared bias and one shared activation to every element, the batched-linear-layer case. It reproduces a loop of gemm_fused calls. gemmkit sizes the single bias vector for one element, not the whole batch.

The slice form: per-element shapes

When the elements differ in shape, or simply do not sit at a fixed stride, use gemm_batched_slice. It takes a slice of BatchProblem, each carrying its own alpha, A, B, beta, and a distinct &mut C view:

#![allow(unused)]
fn main() {
use gemmkit::{gemm_batched_slice, BatchProblem, MatRef, MatMut, Parallelism};

let mut problems: Vec<BatchProblem<'_, f32>> = /* one per product, each its own shape */;
gemm_batched_slice(&mut problems, Parallelism::Rayon(0));
}

Because every C is a distinct &mut, the outputs are pairwise disjoint and cannot alias the inputs by construction. So validation only checks per-element shape agreement and in-bounds strides. Reach for this form when your matrices already live as a Vec of views. Its raw counterpart, gemm_batched_ptr_unchecked over a slice of GemmProblem, takes the same per-element shapes as bare pointers. It serves FFI and adapters that validate their own inputs and may use arbitrary or negative strides. The Unchecked Tier covers both.

How the batch is scheduled

The interesting decision is how the work spreads across cores. The engine makes that decision once per call, from the shared shape and the batch size. There are 3 schedules:

• Batch-parallel. Workers pull disjoint ranges of elements from a shared cursor. Each element runs serially on one worker, cache-hot. This is the whole point of the API. For many small matrices, it pays a single fork/join for the entire batch instead of one per element. It also keeps every core busy on complete GEMMs. This schedule never splits an element across workers, so it is bit-identical to the serial run at any worker count.
• Serial. The whole batch runs on the calling thread, each element single-threaded. The engine picks this when there is too little total work to justify a fork/join.
• Sequential with internal parallelism. For the few-but-large, memory-bound regime, the engine loops the batch and hands each element the full engine parallelism in turn. When an element is big enough to saturate memory bandwidth on its own, spreading it across all cores wins. Running several elements at once instead just thrashes the cache. This schedule runs only for m, n > 1 shapes, whose routes already reduce each output within one worker. So it stays reproducible.

The upshot for you is simple: hand the engine the whole batch and let it choose. Many small independent products are exactly where batching wins over a hand-written loop. A hand-written loop cannot make the fork-join-once-for-all-elements choice. It either parallelizes each tiny GEMM, which is mostly overhead, or runs them serially. For a handful of large products, the batched call converges to what a plain loop already does well. So batching neither helps much nor hurts there.

Determinism spans all 3 schedules. Every element is independent, so the whole batch is reproducible across worker counts. The serial and batch-parallel schedules go further: they are bit-identical across worker counts, because each element runs on exactly one worker. The few-but-large schedule inherits whichever serial-versus-parallel behavior its own per-element route has. A zero-length batch is a no-op.

Like the rest of the API, every batched entry has a _with variant that reuses a caller-owned Workspace, to avoid a per-call allocation. One detail is worth knowing. Under the batch-parallel schedule, the packing cannot go through a single shared Workspace, because concurrent workers would collide on it. So that schedule packs through each worker’s own persistent per-thread pool instead, reused across calls the same way your Workspace is.

Small Shapes and GEMV

The register-blocking driver at the heart of gemmkit needs enough work per output tile. It must amortize packing, cache blocking, and a full MR x NR accumulator.

Some shapes break that premise outright. A matrix-vector product has no tile reuse at all. A k = 4 contraction finishes before the pack pays for itself. An 8 x 8 x 100000 product would spend most of the driver’s effort multiplying padding. For shapes like these, the driver is the wrong tool, so the engine quietly routes around it.

There is nothing to opt into. The same gemm entry (and gemm_i8, gemm_fused, gemm_map) inspects the shape and strides at the top of dispatch. It reroutes to a special-case kernel when one fits the shape, and it falls through to the general driver otherwise. Every reroute stays behind the same public entry.

Each reroute keeps the library’s reproducibility contract: the same call, on the same machine, with the same configuration, returns the same result. Each reroute is also gated by a tuning knob that you can move or disable without recompiling. You never call these paths directly. You benefit from writing the shape naturally, instead of hand-rolling a dot loop.

gemv: the memory-bound edge

A shape with m == 1 or n == 1 is a matrix-vector product. Each output element needs only 2k floating-point operations, but reading it takes k matrix elements. The arithmetic is trivial, so the whole problem is minimizing DRAM traffic. The dedicated gemv path handles both cases with one core routine. It views the matrix as rows x k times a k-vector, transposing the matrix first when m == 1. This path is correct for every layout, and it vectorizes the contiguous ones.

This also covers the degenerate m == n == 1 dot product. A 1 x k view has both strides equal to 1, so it fits the column-major strategy and the row-major strategy at once. The column-major strategy vectorizes over output rows, but there is only one output row here, too few to fill a SIMD register. So the route sends any sweep shorter than one register to the k-vectorizing strategy instead. There is nothing to opt into and nothing to know. A dot product built from MatRef::from_col_major(a, 1, k), the shape a column-major library hands you, runs on the same fast path as the row-major spelling.

2 properties matter to a caller. First, gemv follows the library’s general reproducibility contract. The same call, on the same machine, with the same configuration, always returns the same result. Each output element is reduced in a single pass over k by one worker.

Splitting the rows across workers changes only which worker does the work, not how the work is done. The library does not extend gemv’s guarantee to bitwise agreement across different worker counts. The worker count is part of the configuration, the same as everywhere else in gemmkit.

Second, gemv has its own parallel policy. It is bandwidth-bound, so the worker count comes from a bandwidth model instead of the compute ramp the general driver uses. Past the few cores that saturate DRAM, more workers stop helping, and only add fork/join overhead and shared-cache contention. The count follows a ladder over the bytes a call touches. Below a floor, the matrix fits one core’s private cache, so the call stays serial. Past that floor, the width climbs in steps, and it never reaches the full machine width.

4 knobs expose this policy:

• gemv_parallel_bytes sets the byte floor below which gemv stays single-threaded.
• gemv_tier_step sets how many bytes apart the ladder’s steps sit.
• gemv_thread_cap sets a flat width that replaces the ladder outright.
• gemv_threshold sits alongside the other 3. Since a gemv shape always has min(m, n) == 1, this knob works as an on/off switch rather than a graduated cap.

The small-k path

The contraction k can also be too small for packing to pay off. The threshold is small_k_threshold, with a default of 16 on x86 and 8 on aarch64. At or below that depth, the whole product is a single depth panel, and every packed element would be read exactly once. Packing has nothing to amortize, so the driver’s packing step becomes pure overhead.

The small-k route covers these skinny, low-depth shapes: gevv, rank-k updates, and tall-skinny products. It computes C directly with the family’s microkernel, reading A and B in place, unpacked, in one pass. This route inherits the family’s widen, bias, conjugate, and rounding behavior for free, and it is bit-identical to the serial run for any worker count. It needs a column-major A, meaning unit-stride rows (rsa == 1). When A does not have that layout, packing rarely amortizes at such a small k anyway. The route then defers to the general driver, which still computes the correct result.

The small-m,n path

The mirror case is a tiny output with a long contraction. Both m and n sit far below the microtile, at or below small_mn_dim (default 16 on x86, 32 on aarch64), while k is long. The driver would pad the tiny row and column tiles up to a full microtile. It would then spend most of its work on that padding.

This route instead computes each output as a horizontal dot, C[i,j] = alpha * <A[i,:], B[:,j]> + beta * C[i,j]. It streams SIMD along the contraction with no blocking or orientation machinery. It also register-blocks the small output grid, so several independent FMA chains stay in flight.

The horizontal kernel needs both operands unit-stride along k. That means A’s rows must be contiguous (csa == 1, a row-major A) and B’s columns must be contiguous (rsb == 1, a column-major B). When both hold, the route reads A and B in place with zero copies. This is the fast path.

The 2 most common layouts each miss exactly one side. An all-row-major pair fails rsb, and an all-column-major pair fails csa. When that happens, an internal pre-pack step copies only the failing operand once into k-contiguous scratch, then runs the same horizontal dot over it. That copy reads roughly m*k elements, or n*k for the other operand.

That is a fraction of about 1/n (or 1/m) of the m*n*k work in the product itself. It costs far less than what the horizontal route saves, so a strided small-m,n shape still beats falling back to the driver’s padded microtile.

A flop count understates that copy. The copy does no arithmetic per byte it moves. The dots do about 2. The copy therefore has less to hide memory latency behind, and it takes a much larger share of the time than of the work.

The copy runs across workers on its own, apart from the dots. A small m, n leaves a tiny output grid for the dots to split. The copy splits the depth instead, and a deep contraction makes the depth long. You configure nothing. On the reference machine, a column-major small-m,n shape with a deep contraction got 1.1-3.1x faster (f32, auto width). The larger gains land where the packed operand still fits cache.

The pre-pack tier engages once k clears its own knob, small_mn_pack_min_k (default 16), separate from the zero-copy tier’s small_k_threshold. This path is likewise bit-identical to the serial run at any worker count. It computes each output as one fixed-order reduction on a disjoint tile.

Practical guidance

Prefer the ordinary entries over a hand-written dot loop. Your problem might be a matrix-vector product, a rank-k update, or a grid of small outputs over a long contraction. In each case, gemm already carries a kernel tuned for that shape. It comes with the bandwidth-aware threading and the reproducibility guarantees that a hand-written loop would have to reinvent.

Layout is the one lever in your hands. For the horizontal small-m,n path, a row-major A and a column-major B stream unit-stride and hit the zero-copy fast path. For the small-k path, a column-major A stays on the in-place route. Any other layout still works, but it pays either the small-m,n pre-pack copy or a fall-back to the general driver.

Every threshold named in this chapter is a tunable knob. Each one resolves per call from an argument, a programmatic setter, a GEMMKIT_* environment variable, or a calibrated compile-time default, in that order. If a reroute is miscalibrated for your machine, or if you want to force a shape onto the general driver, see the Tuning Knobs chapter. It covers gemv_threshold, small_k_threshold, small_mn_dim, small_mn_pack_min_k, gemv_parallel_bytes, gemv_tier_step, gemv_thread_cap, and k_stream_max in full. For the special-path internals, meaning why each kernel is shaped the way it is, see the architecture chapter’s Special Paths.

Runtime ISA Dispatch

gemmkit ships one engine. It picks the instruction set it runs on when your program starts, not when you compile it. Take a build made on a laptop and copied to a server. That binary uses the server’s AVX-512 if the server has it. The same binary, run on an older machine, quietly falls back to a narrower kernel. You do not select a backend, gate on a cfg, or rebuild per host. The first GEMM call detects the CPU’s features and caches the winning kernel. Every later call is then a plain indirect call through that cached pointer.

The backend roster

The set of kernels a build carries depends on the target architecture. Which one runs depends on the CPU. From fastest to slowest, the candidates are:

• AVX-512F on x86-64, the widest float kernel. 2 dot-product specializations sit alongside it for the narrow element types: AVX-512 VNNI (vpdpbusd) for i8 -> i32, and AVX-512 BF16 (vdpbf16ps) for bf16. These require the int8 and half features respectively, and the CPU must report the matching feature bit.
• FMA / AVX2 on x86-64, the widen-FMA kernel for machines without AVX-512.
• NEON on aarch64, where SIMD is baseline (every aarch64 CPU has it), so there is nothing to detect at runtime.
• simd128 on wasm32, chosen at compile time rather than runtime (see below).
• scalar, the portable floor. It exists on every target and runs when nothing better is available. A correct, if unaccelerated, result is always reachable.

Tile geometry is the one thing that changes per (element type, ISA) pair. The microkernel computes an MR x NR register tile, sized to the ISA’s vector width. For f32 the shipped tiles are:

All 5 run the same generic float microkernel. Only the tile shape differs. MR is MR_REG * LANES, so a wider vector buys a taller tile. f64 halves the lane count and therefore halves MR (AVX-512F f64 is 16 x 12, and so on). The VNNI and BF16 dot kernels use their own depth-grouped geometry, covered under Element Types. This table is background for reasoning about why a kernel packs and blocks the way it does. It is not a knob you set.

Automatic selection

Each element type owns a single dispatch slot: a OnceLock holding a typed function pointer. On the first call for that type, the selection ladder runs feature detection once and picks the best available kernel. It then stores that kernel’s monomorphized entry points (plain, prepacked, fused), plus the tile geometry. Whichever call happens first pays this one-time cost: the is_x86_feature_detected! probe and the OnceLock initialization. From then on, dispatch is a cached pointer load and an indirect call, with no per-call branching on the ISA. There is no transmute and no atomic pointer juggling behind this, just a typed slot per type.

A consequence worth stating plainly: there is no public API that reports which ISA was selected. The choice is an internal detail of the memoized slot. If you need to be certain a specific kernel is live, do not try to read it back. Pin it instead (next section), and let a mismatch fail loudly.

Pinning a kernel with GEMMKIT_REQUIRE_ISA

Setting the environment variable GEMMKIT_REQUIRE_ISA forces exactly one kernel end to end, instead of auto-selecting. The accepted values (case-insensitive, surrounding whitespace trimmed) are:

The avx512vnni and avx512bf16 pins select the dot kernel for their one narrow type. Everything else runs the plain AVX-512F path, so a mixed workload under one of these pins still runs correctly for its other types.

The contract is panic, not fallback. Dispatch panics, rather than silently running a different kernel, whenever the requested ISA is unavailable. 3 cases count as unavailable:

• the CPU does not report the feature
• the value names an ISA that does not exist on this target architecture (neon on x86, avx512f on aarch64)
• the value is an outright typo

This is deliberate, and it is exactly what you want for CI. A job whose whole purpose is to exercise the AVX-512 VNNI path must not pass by quietly testing the scalar fallback instead. That could happen if a feature flag was misspelled, or if an emulator was misconfigured. gemmkit’s own CI pins each kernel this way: it runs the x86 dot kernels under Intel SDE, NEON on aarch64, and simd128 on wasm. A broken pin then turns into a red build, instead of a false green.

The value is read once, before the first dispatch, and memoized alongside the kernel choice. Set it in the process environment before any GEMM runs. Changing it mid-process has no effect, because the slot is already populated. An unrecognized value is a hard error, precisely so it cannot be mistaken for auto and slip through.

WebAssembly is compile-time

wasm32 has no runtime feature detection, so simd128 is not chosen by probing the machine. It is selected by a compile-time cfg, and the build must actually enable it with -C target-feature=+simd128. Forget the flag, and the wasm build silently uses the scalar floor. Pinning GEMMKIT_REQUIRE_ISA=simd128 turns that silent degradation into an assertion: the build panics if the SIMD path is not live. This is why the wasm CI jobs pin it. See no_std and WebAssembly for the full wasm build story, including the threaded target.

Tuning Knobs

Every heuristic in gemmkit is a named threshold with a shipped default, not a hard-coded constant. This covers questions like when to go parallel, when to pack an operand, and where a shape stops being small. A handful of knobs are split by architecture, where aarch64 needs a different value from the rest. The defaults are good on most hardware. When one is not right for your machine, you can reach it 3 ways without touching the source.

Resolution order

A knob resolves at the point it is read, taking the first of these that is set:

1. Per-call argument. Where a knob has a call-site equivalent, that wins outright. The clearest case is parallelism: the Parallelism argument you pass to gemm overrides any global thread policy. This layer lives in the API, not in tuning.
2. Programmatic setter. gemmkit::tuning::set_*(v) stores a value unconditionally. Once set, later reads never consult the environment again. This is for an application that tunes itself in code. It should win over whatever the deployment environment supplies.
3. Environment variable. GEMMKIT_*. This is the deployment layer. source a profile, for instance one emitted by gemmkit-tune, to retune an already-built binary for a host with no recompile.
4. Compiled default. The calibrated constant, arch-split where it needed to be.

The ordering of setter over env is deliberate. An app that calls the setters has opted out of the environment. An app that wants a deployment profile to apply simply does not call them.

Environment variables are read once, on the first access to that knob, then cached as an atomic. A value set after the first read of a given knob is ignored, so export the profile before the process starts. A GEMMKIT_* var that is set but does not parse as a non-negative integer is treated as a typo, not a silent no-op. gemmkit warns on stderr and falls back to the default. The warning fires once per knob, since the fallback is then cached. It never panics: a perf-knob typo must not crash the process.

The knobs

The table below is the full catalog, across every feature and target configuration. The internal tuning::knob_env_names registry is the source of truth these are drawn from. 2 knobs are feature- or target-gated and only exist when compiled in. Every getter has a matching set_*. The env var name is the getter’s name, upper-cased, with the GEMMKIT_ prefix.

Serial / parallel gate

Pack gates and strides

Special-path thresholds

Scheduler grains

Blocking caps

Deep-contraction and wasm

A note on GEMMKIT_FAST_TEST

You may see GEMMKIT_FAST_TEST in the test harness. It shrinks the correctness sweeps to run faster, and it is a test-suite-only switch. The library itself never reads it, and setting it has no effect on a production GEMM.

Beyond hand-tuning

Setting knobs by hand is for when you already know which one to move. To calibrate the whole set for a specific machine, run the autotuner. It sweeps each knob over a probe-shape set and writes a GEMMKIT_* profile you source before running, with no recompile. That is the subject of the gemmkit-tune chapter.

no_std and WebAssembly

gemmkit’s core does not need an operating system. Turn the default features off, and the crate becomes #![no_std]. It then needs only core and alloc, and depends on nothing else. This makes it usable in kernels, embedded firmware, and WebAssembly. The same code path is also how the wasm SIMD backend gets built. This page covers what a no_std build gives up, what it keeps, and the extra steps a wasm target needs.

The no_std core

The std feature is on by default, as part of default = ["std", "parallel"]. Switch the defaults off, and you get the alloc-only path:

[dependencies]
gemmkit = { version = "0.1", default-features = false }

alloc is always required, because packing scratch is heap-backed in both builds. Beyond that, the crate pulls in nothing at all in this configuration. Each optional feature adds at most one dependency:

• std pulls raw-cpuid, on x86 only, for CPUID cache and feature detection.
• parallel pulls rayon.
• half pulls half.
• complex pulls num-complex.

The int8 and epilogue features add no dependency at all. The element-type features compose freely with no_std. A default-features = false, features = ["half", "int8"] build, for example, is a valid, dependency-free f16/bf16/i8 engine.

Note that parallel implies std, because rayon needs the standard library. So a no_std build is always single-threaded. Everything still compiles and runs. It just runs on one thread.

What changes without std

3 things move from runtime to compile time, or from automatic to explicit:

Feature detection becomes compile-time. With std, x86 dispatch calls is_x86_feature_detected!. It picks the best kernel the running CPU reports. Without std, there is no runtime CPU detection. That capability lives in the std-gated raw-cpuid. So the ISA ladder falls back to cfg!(target_feature = ...) instead, and the build runs whatever its compile-time target features guarantee.

To get an accelerated x86 kernel from a no_std build, you must compile for it. For example, pass -C target-cpu=native or an explicit -C target-feature=+avx512f. Otherwise you get the scalar floor. On aarch64 and on wasm, this is already how selection works, so nothing is lost there.

The env knobs are off. Reading an environment variable needs std. Without it, GEMMKIT_REQUIRE_ISA is never consulted (dispatch always auto-selects), and every GEMMKIT_* tuning knob resolves straight to its compiled default. The programmatic tuning::set_* setters still work, so you retune a no_std build in code rather than through the environment. See Tuning Knobs for the setter layer.

A per-call workspace replaces the pool. The default thread-local packing pool is a std construct. Without std, there is no pool. Each call instead allocates a fresh Workspace for its scratch, and frees it on return. That is correct, but it allocates on every call. To reach a zero-allocation steady state, create a Workspace once, and thread it through the *_with entries: gemm_with, and the _with variant of every family. After the first sufficiently large call, those entries reuse the buffer with no further heap traffic.

Building for WebAssembly

wasm32 has no runtime feature detection, so the simd128 backend is selected by a compile-time cfg. The build must enable that target feature explicitly. If you forget it, the wasm build still compiles and runs correctly, but only on the scalar floor, which is much slower. Pass the flag through RUSTFLAGS:

RUSTFLAGS="-C target-feature=+simd128" \
  cargo build --target wasm32-wasip1 --no-default-features --features std

To run the result, you need a wasm runtime. gemmkit’s CI uses wasmtime, so point Cargo’s target runner at it. Sometimes you want to be certain the SIMD path is actually live, rather than silently falling back to scalar. Pin the ISA for that: GEMMKIT_REQUIRE_ISA=simd128 turns a missing +simd128 into a panic, instead of a quiet fallback. That is exactly what a test job wants. This pin needs std, which is already on in the wasm builds.

RUSTFLAGS="-C target-feature=+simd128" \
CARGO_TARGET_WASM32_WASIP1_RUNNER="wasmtime --env GEMMKIT_REQUIRE_ISA=simd128" \
  cargo test --target wasm32-wasip1 --no-default-features --features std

Baseline wasm32-wasip1 has no threads. If you build with parallel on for a baseline wasm target, gemmkit does not trap. An internal guard makes rayon unusable there, so Parallelism::Rayon(_) degrades to the serial loop instead. A portable wasm binary can therefore carry the parallel feature and simply run single-threaded, with no target-specific build.

Threaded wasm

Real multithreading on wasm needs the threads-capable target and the matching feature:

RUSTFLAGS="-C target-feature=+simd128" \
CARGO_TARGET_WASM32_WASIP1_THREADS_RUNNER="wasmtime -W threads=y -W shared-memory=y -S threads=y" \
  cargo test --target wasm32-wasip1-threads \
  --no-default-features --features std,parallel,wasm_threads

The wasm_threads feature implies parallel, and it targets wasm32-wasip1-threads. It turns on gemmkit’s dedicated wasm rayon pool. Because a wasm runtime cannot report a core count, the pool’s width is not auto-derived. It comes from the GEMMKIT_WASM_THREADS knob instead, default 8, which both caps the auto worker count and sizes the pool. Set it to match the number of workers your runtime actually provisions. Everything else behaves exactly as on a native threaded build: blocking, the job list, and reproducibility.

The Unchecked Tier

The safe entries (gemm, gemm_fused, and the rest) validate their inputs before touching memory. Shapes must agree. Every strided view must stay inside its slice. The output must address each element once, and it must not overlap the inputs.

Under every one of those checks sits the same engine. It is reached through a raw-pointer, isize-stride interface with no checks at all. That is the unchecked tier. It exists for callers who already hold the invariants that the safe API would otherwise re-derive.

Who it is for

3 kinds of caller live here. Adapters over other matrix libraries, such as ndarray, nalgebra, and faer, already have a validated pointer and strides straight out of the host type. Re-checking bounds would be redundant work on data the library already guarantees. FFI callers arriving from C or another language have a pointer and strides, and no Rust slice to bound-check against. Custom matrix types that a codebase owns can lower to pointers and call the engine directly, instead of copying into a MatRef. In each case, the caller is the party that knows the memory is valid, so the checks move to where that knowledge lives.

If none of that describes you, use the safe API. The unchecked tier is not faster for a single call. The validation cost is cheap relative to the multiply itself. It exists to let a caller that already owns the invariants avoid proving them again.

The catalog

Every safe entry has a raw twin, named by appending _unchecked, and most also offer a _with form that takes a caller-owned workspace (next section). The full raw surface, by family:

The pointer-array batched form is worth calling out. gemm_batched_ptr_unchecked takes a slice of GemmProblem<T>, each with its own shape and its own pointers. So a batch can mix sizes, and it can scatter its operands anywhere in memory. It has no safe counterpart of the same shape. Expressing “an array of independent raw problems” is precisely what the raw tier is for. The nalgebra and faer adapters build their batched GEMM on top of it.

The safety contract

Calling into the unchecked tier means signing, per call, for what the safe API would otherwise check:

• Valid pointers and strides. For every (i, j) implied by the dimensions and strides, a and b are valid for reads and c is valid for read and write. Nothing bounds-checks this. An out-of-range stride is undefined behavior, not a panic.
• A uniquely-addressed output. C’s strides must map every distinct (i, j) to a distinct location. The parallel driver assumes output tiles are disjoint, and it writes them concurrently. A self-aliasing C (for example rsc == 0) would then be a data race. The inputs may alias themselves freely, since they are only read. So a broadcast, zero-stride, A or B is fine.
• No overlap between C and A/B. The output is written. If it overlapped an input, the result would be garbage.

One relaxation comes with the territory. When beta == 0, the output is not read, so C need not be initialized. One capability the safe API withholds is available here too: negative strides, and pointers into the middle of a buffer, are both allowed. A reversed view (rs < 0), or an operand addressed backward from its last element, is exactly the kind of layout the safe MatRef refuses. The raw engine accepts it instead. This is why adapters over libraries that produce reversed strides forward to this tier.

Reusing a workspace

Each raw entry comes in 2 allocation flavors. The plain form, gemm_unchecked, borrows the thread-local packing pool. It allocates at most once per thread. The _with form, gemm_unchecked_with, takes a &mut Workspace that you own instead:

#![allow(unused)]
fn main() {
use gemmkit::{Workspace, Parallelism};

let mut ws = Workspace::new();
// each iteration reuses `ws`; after the first large call it does no heap work
for _ in 0..iters {
    // SAFETY: pointers/strides valid, c uniquely addressed, c disjoint from a/b
    unsafe {
        gemmkit::gemm_unchecked_with(
            &mut ws, m, k, n,
            1.0_f32, a, rsa, csa, b, rsb, csb, 0.0_f32, c, rsc, csc,
            Parallelism::Serial,
        );
    }
}
}

The workspace grows to fit the largest problem it has served. It reuses that allocation thereafter, so a hot loop of GEMMs reaches zero steady-state allocation. This is the mechanism that no_std builds rely on for reuse, since they have no thread-local pool. It is equally useful under std, for real-time or latency-sensitive loops where you want allocation off the hot path.

A worked example: a custom tile type

Suppose your code already carries its own dense row-major matrix, and you want to multiply two of them without copying into a MatRef:

#![allow(unused)]
fn main() {
use gemmkit::{gemm_unchecked, Parallelism};

// a dense row-major matrix the caller already owns
struct Tile {
    data: Vec<f32>,
    rows: usize,
    cols: usize,
}

// c = a * b for row-major tiles
fn matmul(a: &Tile, b: &Tile, c: &mut Tile) {
    assert_eq!(a.cols, b.rows);
    assert_eq!(a.rows, c.rows);
    assert_eq!(b.cols, c.cols);
    // row-major: row stride = cols, column stride = 1
    // SAFETY: shapes checked above; each tile owns a dense rows*cols buffer, so
    // every addressed element is in bounds; c is a distinct &mut, so it cannot
    // alias a or b, and a dense layout addresses each (i, j) once
    unsafe {
        gemm_unchecked(
            a.rows, a.cols, b.cols,
            1.0_f32,
            a.data.as_ptr(), a.cols as isize, 1,
            b.data.as_ptr(), b.cols as isize, 1,
            0.0_f32,
            c.data.as_mut_ptr(), c.cols as isize, 1,
            Parallelism::Serial,
        );
    }
}
}

The assert_eq! shape checks and the &mut Tile borrow together discharge the whole contract. Dense storage makes every offset land in bounds, and it makes every (i, j) distinct. The exclusive borrow of c then rules out overlap with a or b. This is the pattern to reach for: prove the invariants at the boundary of your own type, then hand raw pointers to the engine.

The adapters are the reference

The cleanest examples of doing this well are the adapter crates themselves. Each one pulls a pointer and strides out of a native view: C-order, F-order, general strides, or reversed strides, all with no copies. Each then forwards to the *_unchecked engine, with a short safety argument at each call site. If you are wrapping a matrix type of your own, read one of the adapter chapters and mirror its structure. The nalgebra chapter is a good place to start. For the prepacked entries in the catalog above, the fixed-weight reuse pattern they serve is covered in Prepacked Operands.

Using gemmkit with ndarray

gemmkit-ndarray is a thin bridge between ndarray’s two-dimensional arrays and the gemmkit engine. It does no numerical work of its own. Each entry takes an ArrayBase. It reads the base pointer and the 2 axis strides straight out of it, and hands those raw parts to gemmkit’s unchecked engine. The whole crate is a stride-plumbing layer. Everything gemmkit knows how to do, such as runtime ISA selection, cache blocking, and reproducible parallelism, applies unchanged. The arrays never get reshaped or copied on the way in.

The entries accept &ArrayBase<S, Ix2> for any storage S: Data. Both an owned &Array2<T> and a borrowed ArrayView2<T> work, along with ArcArray, CowArray, and slices of any of them. The only internal helper is worth seeing. It is the whole of the adapter’s data extraction:

gemmkit-ndarray/src/common.rs:

#![allow(unused)]
fn main() {
pub(crate) fn dims_strides<T, S: Data<Elem = T>>(
    a: &ArrayBase<S, Ix2>,
) -> (usize, usize, isize, isize) {
    let (r, c) = a.dim();
    let s = a.strides();
    (r, c, s[0], s[1])
}
}

That (rows, cols, row_stride, col_stride) tuple, plus a.as_ptr(), is everything gemmkit needs. The strides are signed isize, so a negative (reversed) stride forwards just like a positive one.

Adding it to a project

2 crates, not 3. The adapter re-exports everything its own signatures name, so a direct gemmkit dependency is not part of the normal setup:

[dependencies]
gemmkit-ndarray = "0.1"
ndarray = "0.17.1"

gemmkit_ndarray re-exports everything a caller needs, so a direct gemmkit dependency is rarely necessary:

• The Parallelism selector and the Workspace type every _with variant takes.
• The fused selectors Bias and Activation.
• The prepacked handles PackedLhs and PackedRhs.
• The requantization parameters Requantize and RequantScale.
• The element-type bounds GemmScalar, FusedScalar, MapScalar, and ComplexScalar. Name these when you write a wrapper generic over an entry.
• The element types f16, bf16, Complex, c32, and c64, each under its own feature. This keeps half and num-complex out of your manifest too.
• The tuning module.

Reach for tuning through the adapter. Do not add a separate gemmkit dependency of your own for it. The tuning knobs are process-global atomics. A second, separately resolved gemmkit would give you a set of atomics the adapter never reads.

Every feature on gemmkit-ndarray forwards directly to the same-named feature on gemmkit. Turn on a capability here, and the matching entry points light up:

• parallel (default): rayon multithreading.
• wasm_threads: threading on wasm32-wasip1-threads. Implies parallel.
• half: f16 / bf16 inputs with f32 accumulation.
• complex: Complex<f32> / Complex<f64> matrices.
• int8: i8 inputs accumulating into i32.
• epilogue: fused bias / activation, i8 / u8 requantization, and a user per-element map.

The default is ["parallel"]. The advanced page covers the feature-gated families, such as gemm_cplx, gemm_i8, and gemm_fused. The minimum ndarray version is 0.17.1.

The core entries

3 functions cover the plain real path. dot is the convenience entry. It multiplies A * B into a freshly allocated row-major Array2, the way ndarray’s own .dot() reads.

#![allow(unused)]
fn main() {
use ndarray::array;

let a = array![[1.0_f32, 2.0], [3.0, 4.0]];
let b = array![[5.0_f32, 6.0], [7.0, 8.0]];
let c = gemmkit_ndarray::dot(&a, &b);
assert_eq!(c, array![[19.0, 22.0], [43.0, 50.0]]);
}

dot is generic over T: GemmScalar, which is f32 and f64 unconditionally, plus f16 and bf16 when the half feature is on. It parallelizes with Parallelism::default() and allocates its own output. Use it for a one-off product where you do not already own the destination.

gemm writes the general form C <- alpha*A*B + beta*C in place. This is where alpha, beta, an existing accumulator, and an explicit Parallelism come in. Its signature is:

#![allow(unused)]
fn main() {
pub fn gemm<T, S1, S2, SC>(
    alpha: T,
    a: &ArrayBase<S1, Ix2>,
    b: &ArrayBase<S2, Ix2>,
    beta: T,
    c: &mut ArrayBase<SC, Ix2>,
    par: Parallelism,
)
where
    T: GemmScalar,
    S1: Data<Elem = T>,
    S2: Data<Elem = T>,
    SC: DataMut<Elem = T>;
}

The output binds SC: DataMut, so C is a &mut Array2 or an ArrayViewMut2 and, like the inputs, may carry any layout. Here A is a row-major buffer transposed into a column-major view with no copy. The multiply runs single-threaded:

#![allow(unused)]
fn main() {
use gemmkit_ndarray::Parallelism;
use ndarray::{Array2, array};

// row-major storage, transposed into a column-major view with no copy
let a = Array2::from_shape_vec((2, 2), vec![1.0_f32, 2.0, 3.0, 4.0])
    .unwrap()
    .reversed_axes();
let b = Array2::from_elem((2, 2), 1.0_f32);
let mut c = Array2::zeros((2, 2));
gemmkit_ndarray::gemm(1.0, &a, &b, 0.0, &mut c, Parallelism::Serial);
assert_eq!(c, array![[4.0, 4.0], [6.0, 6.0]]);
}

Layouts that cost nothing

Because the adapter only reads strides, any two-dimensional view that ndarray can express forwards without a copy. That includes:

• the standard C-order (row-major) layout
• an F-order (column-major) view from .reversed_axes(), .t(), or an array built with .f()
• a windowed .slice(...) view with non-unit strides
• a reversed view from a negative-step slice such as s![..;-1, ..], which produces a negative row stride

The destination C is just as free. Array2::zeros((m, n).f()) gives a column-major output, and gemm fills it directly.

“Zero-copy” here means the adapter never copies to normalize a layout. gemmkit’s engine still packs operands into its own scratch buffers when the microkernel needs contiguous panels. That internal packing is part of the algorithm, not a materialization of a transposed input. The point is that you never pay a to_owned() or a manual transpose to satisfy the call, whatever your arrays look like.

Panics: shapes, not aliasing

The adapter validates shapes and nothing else. Each entry asserts that the inner dimensions line up and that C matches the product. On a mismatch it panics with a gemmkit-ndarray: message that names the offending dimensions, for instance A.cols (k) != B.rows (kb), A.rows (m) != C.rows (cm), or B.cols (n) != C.cols (cn). A dimension mismatch is the only reason a plain gemm or dot panics.

The adapter does not check aliasing at runtime, and it does not need to. C arrives as &mut ArrayBase<SC, _>, an exclusive borrow. The type system already guarantees it cannot overlap the shared & borrows of A and B. That is exactly the precondition gemmkit’s _unchecked engine asks its caller to uphold, and the &mut signature upholds it for free. The fused entries add 1 more runtime check, for a bias slice overlapping C. See the advanced page for that check.

Choosing parallelism

Parallelism is re-exported by the adapter. Parallelism::Serial runs on the calling thread. Parallelism::Rayon(n) uses a rayon pool of at most n threads, and Rayon(0) auto-detects the machine’s core count. Parallelism::default() is Rayon(0), which is what dot uses, so dot is parallel out of the box. The threaded paths need the parallel feature, which is on by default. With that feature off, treat every call as serial.

Blocking and job order do not depend on the thread count, so a fixed input and configuration give a reproducible result. Serial and parallel runs agree bit-for-bit today, because both walk the same blocking and the same kernel. That agreement is a property of the current implementation, not a promise across every configuration. The reasoning behind the thread counts lives in Parallelism in Practice.

Reusing a workspace

Every allocating entry borrows scratch from gemmkit’s internal thread-local pool for the duration of the call. A lone gemm call never leaks an allocation into your steady state. When you run a hot loop of similar products, the _with variants let you own that scratch instead. Pass a &mut Workspace as the first argument. It grows once to the largest size the loop needs, then gets reused with no further allocation.

#![allow(unused)]
fn main() {
use gemmkit_ndarray::{Parallelism, Workspace};
use ndarray::Array2;

let mut ws = Workspace::new();
let par = Parallelism::default();
for &(m, k, n) in &[(256, 256, 256), (512, 128, 512)] {
    let a = Array2::<f32>::zeros((m, k));
    let b = Array2::<f32>::zeros((k, n));
    let mut c = Array2::<f32>::zeros((m, n));
    // reuses ws across iterations; allocates at most once
    gemmkit_ndarray::gemm_with(&mut ws, 1.0, &a, &b, 0.0, &mut c, par);
}
}

gemm_with is identical to gemm apart from the leading workspace, and it produces the same result. Every family in the adapter has a matching _with twin. The pattern carries over to the integer, complex, fused, batched, and prepacked entries covered next. If you multiply a fixed weight matrix against a stream of activations, the workspace pairs naturally with the prepacked-operand path on the advanced page.

ndarray Adapter Advanced Usage

Beyond the plain real product, the adapter mirrors the whole of gemmkit’s surface:

• integer GEMM
• requantized quantized-inference output
• complex products with optional conjugation
• fused bias and activation
• a user-supplied per-element map
• batched multiplication over rank-3 arrays
• prepacked operands

The adapter gates each family behind the Cargo feature named for it. Each family also keeps the shape of the plain entries from the getting-started page. Each reads strides straight from the arrays, forwards to gemmkit, and panics only on a dimension mismatch. The fused entries add 1 more panic case, for a bias slice that overlaps C. Every entry also has a _with twin that threads a caller-owned Workspace.

Integer GEMM (int8)

gemm_i8 multiplies i8 inputs into an i32 accumulator: C(i32) <- alpha*A(i8)*B(i8) + beta*C, with alpha, beta, and C all i32. It is a separate entry from gemm because the input and output element types differ. Arithmetic wraps on overflow. This is the conventional integer-GEMM contract. dot_i8 is the convenience twin. It returns a fresh Array2<i32>.

#![allow(unused)]
fn main() {
use gemmkit_ndarray::{Parallelism, dot_i8, gemm_i8};
use ndarray::Array2;

let a = Array2::<i8>::zeros((16, 12));
let b = Array2::<i8>::zeros((12, 10));

// i8 inputs, i32 accumulator
let c: Array2<i32> = dot_i8(&a, &b);

// general form with i32 alpha/beta into an existing accumulator
let mut acc = Array2::<i32>::zeros((16, 10));
gemm_i8(2, &a, &b, 1, &mut acc, Parallelism::Serial);
}

Requantized output (int8 + epilogue)

Quantized inference rarely wants the raw i32 accumulator. It wants an 8-bit tensor back. gemm_i8_requant fuses the multiply and the requantize into 1 pass. It folds the i32 accumulator to an i8 output without ever materializing the full m*n intermediate. There is no alpha parameter, because it folds into the scale, and no beta parameter, because accumulating into a quantized output is ill-defined. The parameters live in a Requantize:

#![allow(unused)]
fn main() {
use gemmkit_ndarray::{Parallelism, RequantScale, Requantize, gemm_i8_requant, gemm_i8_requant_u8};
use ndarray::Array2;

let a = Array2::<i8>::zeros((16, 12));
let b = Array2::<i8>::zeros((12, 10));

// i8 output in [-128, 127], one per-tensor scale, per-row bias (length A.rows)
let bias: Vec<i32> = vec![0; 16];
let mut c = Array2::<i8>::zeros((16, 10));
let req = Requantize {
    scale: RequantScale::PerTensor(0.05),
    zero_point: -7,
    bias: Some(&bias),
};
gemm_i8_requant(&a, &b, req, &mut c, Parallelism::default());

// u8 output in [0, 255], per-channel scales, no bias
let scales: Vec<f32> = vec![0.02; 16]; // one per output row / channel
let mut cu = Array2::<u8>::zeros((16, 10));
gemm_i8_requant_u8(
    &a,
    &b,
    Requantize { scale: RequantScale::PerRow(&scales), zero_point: 128, bias: None },
    &mut cu,
    Parallelism::default(),
);
}

The output is clamp(zero_point + round_ne(scale * (accumulator + bias[i])), LO, HI) with round-half-to-even, where scale is the per-tensor value or the per-row scale_i. The u8 variant is the ONNX-QLinearMatMul-style activation: identical to gemm_i8_requant apart from the output domain [0, 255] and the zero_point band. Both entries reject:

• a non-finite or non-positive scale, per-tensor or per-row
• a per-row scale or bias whose length is not A.rows
• a slice that overlaps C
• a zero_point outside the entry’s domain

Complex GEMM (complex)

Complex products get their own entries, because the 2 conjugation flags do not fit the homogeneous real signature. gemm_cplx computes C <- alpha*op(A)*op(B) + beta*C, over Complex<f32> or Complex<f64>. op(A) is conj(A) when conj_a is set. op(B) is conj(B) when conj_b is set. dot_cplx is the non-conjugated convenience.

#![allow(unused)]
fn main() {
use gemmkit_ndarray::{Complex, Parallelism, dot_cplx, gemm_cplx};
use ndarray::Array2;

type C = Complex<f64>;
let a = Array2::<C>::from_elem((8, 6), Complex::new(0.0, 0.0));
let b = Array2::<C>::from_elem((6, 5), Complex::new(0.0, 0.0));

// plain A*B
let c = dot_cplx(&a, &b);

// conjugate A, accumulate into an existing C
let mut acc = Array2::<C>::from_elem((8, 5), Complex::new(0.0, 0.0));
gemm_cplx(
    Complex::new(1.0, 0.0),
    &a,
    true,  // conj_a
    &b,
    false, // conj_b
    Complex::new(0.0, 0.0),
    &mut acc,
    Parallelism::Serial,
);
}

With complex and epilogue both on, gemm_cplx_fused adds an optional Bias (added verbatim, never conjugated) in the same pass. It takes no activation parameter. An ordering activation like ReLU is undefined on complex numbers.

Fused bias, activation, and maps (epilogue)

gemm_fused computes C <- act(alpha*A*B + beta*C + bias) in 1 pass. The bias is an optional Bias::PerRow (length A.rows) or Bias::PerCol (length B.cols). The activation is an optional Relu or LeakyRelu(slope), applied last. With both set to None, gemm_fused is exactly gemm.

#![allow(unused)]
fn main() {
use gemmkit_ndarray::{Activation, Bias, Parallelism, gemm_fused, gemm_map};
use ndarray::Array2;

let a = Array2::<f32>::zeros((12, 9));
let b = Array2::<f32>::zeros((9, 7));

// C <- ReLU(A*B + bias) in one pass; PerRow bias has length A.rows
let bias: Vec<f32> = vec![0.0; 12];
let mut c = Array2::<f32>::zeros((12, 7));
gemm_fused(
    1.0, &a, &b, 0.0, &mut c,
    Some(Bias::PerRow(&bias)),
    Some(Activation::Relu),
    Parallelism::default(),
);

// arbitrary per-element closure f(value, row, col); here a relu6
let f = |v: f32, _r: usize, _c: usize| v.max(0.0).min(6.0);
let mut c2 = Array2::<f32>::zeros((12, 7));
gemm_map(1.0, &a, &b, 0.0, &mut c2, &f, Parallelism::default());
}

Bias and Activation are re-exported from gemmkit_ndarray, so you need not name gemmkit for them. For f32/f64, gemm_fused is bit-identical to gemm followed by the same scalar map, for every shape. For f16/bf16, the epilogue runs in f32 before the single narrowing. This is more precise than a separate narrow map, so the result is not bitwise-equal to gemm-then-map for those types.

gemm_map is the general per-element extension point. The closure f(value, row, col) sees each output element at its final value, with (row, col) in the user frame of C. gemmkit calls it exactly once per element. It costs 1 indirect call per element. Prefer gemm_fused for a plain bias or activation, because it vectorizes. Reach for gemm_map instead for GELU, sigmoid, clamps, or position-dependent transforms. Here, T is f32/f64 only.

Batched GEMM

This is the only operation with no plain-gemm analogue, and no counterpart in the sibling adapters. It is a stack of independent products on a rank-3 Array3, with the batch on axis 0. a is (batch, m, k), b is (batch, k, n), and c is (batch, m, n). Axis 0 is each operand’s batch stride. Axes 1 and 2 are the element strides. gemm_batched parallelizes across the batch. Each element runs on 1 worker, so the result reproduces a loop of gemm calls exactly.

#![allow(unused)]
fn main() {
use gemmkit_ndarray::{Parallelism, dot_batched, gemm_batched};
use ndarray::Array3;

let a = Array3::<f32>::zeros((32, 8, 5)); // (batch, m, k)
let b = Array3::<f32>::zeros((32, 5, 6)); // (batch, k, n)

// stack of products
let c = dot_batched(&a, &b); // (32, 8, 6)

// general form into an existing accumulator
let mut acc = Array3::<f32>::zeros((32, 8, 6));
gemm_batched(0.7, &a, &b, 1.3, &mut acc, Parallelism::default());
}

The adapter reads only strides, so a permuted-axes or otherwise general-stride 3-D view forwards without a copy. For example, a.view().permuted_axes([0, 2, 1]) turns a (batch, k, m) buffer into a (batch, m, k) view that batches straight through.

Under epilogue, gemm_batched_fused applies 1 shared Bias/Activation to every element of the stack. This is the batched-linear-layer case. The bias is sized for a single element (PerRow length m, PerCol length n), not the whole batch.

Prepacked operands

When one operand is fixed and the other streams, pack the fixed side once and reuse it. This skips the per-call repack. prepack_rhs returns a PackedRhs<T> for a reused B, consumed by gemm_packed_b. prepack_lhs returns a PackedLhs<T> for a reused A, consumed by gemm_packed_a. The prepack functions read strides directly, so B or A may have any layout.

Each has 1 orientation constraint. gemm_packed_b needs a column-major-ish C (|col stride| >= |row stride|). gemm_packed_a needs a row-major-ish C (|col stride| <= |row stride|). The other orientation would swap the operands and invalidate the packed handle, which gemmkit rejects. Use plain gemm for the layout that does not fit.

The fused twins gemm_packed_b_fused and gemm_packed_a_fused accept the same handles plus a bias and activation. This is exactly the fixed-weight inference layer. The example below packs a weight matrix once as the LHS, and reuses it across inference steps. It folds in a per-output-channel bias and a ReLU:

#![allow(unused)]
fn main() {
use gemmkit_ndarray::{Activation, Bias, Parallelism, gemm_packed_a_fused, prepack_lhs};
use ndarray::Array2;

let (out, in_features) = (256usize, 512usize);

// pack the fixed weight W: (out, in) once
let w = Array2::<f32>::zeros((out, in_features));
let packed = prepack_lhs(&w);
let bias: Vec<f32> = vec![0.0; out]; // per-output-channel, length C.rows

// each inference step: activations x (in, batch) -> y (out, batch)
let batch = 32;
let x = Array2::<f32>::zeros((in_features, batch));
let mut y = Array2::<f32>::zeros((out, batch)); // row-major (packed_a orientation)
gemm_packed_a_fused(
    1.0,
    &packed,
    &x,
    0.0,
    &mut y,
    Some(Bias::PerRow(&bias)),
    Some(Activation::Relu),
    Parallelism::default(),
);
}

Pair the packed handle with the _with workspace variant, gemm_packed_a_fused_with. A steady inference loop then allocates nothing after the first call. You specify the bias axis in the user frame. The packed path forwards it unflipped for gemm_packed_b_fused, and lets the core flip it for gemm_packed_a_fused. PerRow always means “1 value per output row,” regardless of which operand you packed.

This adapter versus ndarray’s own product

ndarray already multiplies matrices: .dot() for the plain product and general_mat_mul for the in-place alpha/beta form. For a one-off f32/f64 product with no extra requirements, those functions work well and pull in 1 less dependency. There is no reason to route through gemmkit out of habit.

Reach for this adapter when you want what ndarray’s built-in path does not offer. gemmkit picks the fastest instruction set on the machine it runs on, at runtime. It does not bake 1 choice in at compile time (see Runtime ISA Dispatch). It brings the wider surface this page covers: fused bias and activation, i8 and requantized inference, complex with conjugation, batched products, and prepacking. It also exposes tuning knobs, plus an install-time autotuner that calibrates blocking to the deployment machine (see Tuning Knobs).

None of that changes the arrays you pass or the results you get back. The adapter is the same zero-copy stride plumbing throughout. It only widens what you can ask for.

Using gemmkit with nalgebra

gemmkit-nalgebra lets you drive the gemmkit engine straight from nalgebra matrices without copying them first. It targets nalgebra 0.35 and accepts &Matrix<T, R, C, S> for any storage S: RawStorage<T, R, C>. Owned DMatrix, static SMatrix, and every view or slice type all qualify. The adapter reads the matrix’s data pointer and its 2 strides. It hands them to gemmkit’s raw engine. The engine does not reshape the input, transpose it into a scratch buffer, or duplicate it on the way in.

nalgebra’s natural layout is column-major, which is also gemmkit’s preferred orientation, so the common case is the fast case. Row-major and general-stride views work too, with the same zero-copy behavior, because the engine reads strides directly instead of assuming one layout.

Adding it to a project

2 crates go into Cargo.toml. The adapter re-exports everything its own signatures name, so a direct gemmkit dependency is not part of the normal setup.

[dependencies]
gemmkit-nalgebra = "0.1"
nalgebra = "0.35"

gemmkit_nalgebra re-exports several types, so a project does not need a direct gemmkit dependency for them:

• the Parallelism selector and the Workspace every _with variant takes
• the fused selectors Bias and Activation
• the prepacked handles PackedLhs and PackedRhs
• the requantization parameters Requantize and RequantScale
• the element-type bounds GemmScalar, FusedScalar, MapScalar, and ComplexScalar. Name these when you write a wrapper generic over an entry
• the element types f16, bf16, Complex, c32, and c64, gated behind their features, so half and num-complex also stay out of your manifest
• the tuning module

Reach for tuning through the adapter instead of through your own gemmkit dependency. The tuning knobs are process-global atomics. A second, separately resolved gemmkit crate would give you a set of atomics the adapter never reads.

The default feature set enables parallel, which turns on rayon-based threading in the engine (gemmkit/parallel). Every adapter feature is a thin forward to the same-named feature on gemmkit:

• half adds f16/bf16 inputs
• complex adds Complex<f32>/Complex<f64>
• int8 adds the i8 -> i32 path
• epilogue adds fused bias/activation and the per-element map
• wasm_threads layers parallel onto wasm32-wasip1-threads

The feature-gated entries are covered in nalgebra Adapter Advanced Usage. This page stays on the always-available real-scalar surface.

The 3 real-scalar entries

The base surface is 3 functions, all generic over GemmScalar (which is f32 and f64 always, plus f16 and bf16 under the half feature). gemm is the accumulating multiply, gemm_with is the same call reusing a caller-owned workspace, and dot is a convenience wrapper that allocates its result.

Here is gemm verbatim, from gemmkit-nalgebra/src/float.rs:

#![allow(unused)]
fn main() {
pub fn gemm<T, R1, C1, S1, R2, C2, S2, RC, CC, SC>(
    alpha: T,
    a: &Matrix<T, R1, C1, S1>,
    b: &Matrix<T, R2, C2, S2>,
    beta: T,
    c: &mut Matrix<T, RC, CC, SC>,
    par: Parallelism,
) where
    T: GemmScalar,
    R1: Dim,
    C1: Dim,
    S1: RawStorage<T, R1, C1>,
    R2: Dim,
    C2: Dim,
    S2: RawStorage<T, R2, C2>,
    RC: Dim,
    CC: Dim,
    SC: RawStorageMut<T, RC, CC>,
{
    gemm_common(None, alpha, a, b, beta, c, par);
}
}

The 10 generic parameters look heavy, but they say one simple thing. A, B, and C may each be any nalgebra matrix or view, with independent row-dimension, column-dimension, and storage types. A and B are read through RawStorage. C needs RawStorageMut, because gemm writes into it in place. The operation is C <- alpha*A*B + beta*C.

gemm_with takes the same arguments, preceded by a &mut Workspace. dot(a, b) -> DMatrix<T> computes A*B into a freshly allocated column-major matrix. It calls gemm internally with beta == 0, so the fresh buffer is never read before it is overwritten.

A first multiply with DMatrix

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::Parallelism;
use nalgebra::DMatrix;

let a = DMatrix::from_row_slice(2, 2, &[1.0_f32, 2.0, 3.0, 4.0]);
let b = DMatrix::from_row_slice(2, 2, &[5.0_f32, 6.0, 7.0, 8.0]);

// dot: A*B into a fresh column-major DMatrix
let c = gemmkit_nalgebra::dot(&a, &b);
assert_eq!(c, DMatrix::from_row_slice(2, 2, &[19.0, 22.0, 43.0, 50.0]));

// gemm: accumulate C <- alpha*A*B + beta*C in place
let mut acc = DMatrix::<f32>::zeros(2, 2);
gemmkit_nalgebra::gemm(1.0, &a, &b, 0.0, &mut acc, Parallelism::default());
assert_eq!(acc, c);
}

dot is the right tool when you just want the product and are happy with a DMatrix back. gemm is the tool when you already own the destination. Use it to scale the destination (beta), scale the product (alpha), or avoid the allocation dot performs. dot always returns a DMatrix<T>, even for static inputs, because the wrapper only knows the output dimensions at the value level.

Static and mixed-shape matrices

Static matrices go through the same functions with no special casing. Because the row, column, and storage generics are independent per operand, a static A can multiply a dynamic B. gemm can write into a static &mut SMatrix output just as readily as into a DMatrix.

#![allow(unused)]
fn main() {
use nalgebra::{DMatrix, Matrix2, SMatrix};

// static x static
let a = Matrix2::new(1.0_f32, 2.0, 3.0, 4.0);
let b = Matrix2::new(5.0_f32, 6.0, 7.0, 8.0);
let c = gemmkit_nalgebra::dot(&a, &b); // -> DMatrix<f32>
assert_eq!(c[(0, 0)], 19.0);

// static A x dynamic B: the independent Dim generics allow it
let a34 = SMatrix::<f64, 3, 4>::from_fn(|i, j| (i as f64) - 0.5 * (j as f64) + 1.0);
let b = DMatrix::<f64>::from_element(4, 2, 0.25);
let c = gemmkit_nalgebra::dot(&a34, &b); // -> DMatrix<f64>, shape 3x2
}

Layouts and zero-copy

The adapter never copies an operand. It pulls (rows, cols, row-stride, col-stride) out of the matrix and forwards the pointer plus strides to the engine. The source layout only decides which stride pair the engine sees. It never decides whether an allocation happens. A column-major DMatrix is read in place, and so is a row-major slice built with from_slice_with_strides. A non-contiguous stepped view, such as every other row of a larger matrix, is read in place too. nalgebra reports strides as non-negative element counts, and the adapter widens them to the signed strides the engine expects.

Whatever internal packing the engine does to feed its microkernel is independent of the source layout and happens regardless of where the data came from. That is a property of gemmkit, not a copy the adapter introduces. If you want to eliminate the repeated internal repacking of a reused operand, use the prepacked-operand path instead. It is described in the advanced page.

When it panics

The entries validate shapes and panic on a mismatch, with the offending dimensions in the message. gemm (and its siblings) check 3 equalities before touching any memory: A.cols == B.rows, A.rows == C.rows, and B.cols == C.cols. A mismatch on the inner dimension, for instance, aborts with gemmkit-nalgebra: A.cols (k) != B.rows (kb) instead of reading out of bounds. The output matrix must therefore already have the right shape. gemm writes into it but does not resize it. dot, which allocates the output itself, can only fail on the inner-dimension check.

Choosing parallelism

Every call takes a Parallelism as its last argument. There are 2 variants. Parallelism::Serial runs single-threaded, and Parallelism::Rayon(n) runs on rayon with at most n threads, where Rayon(0) auto-detects the thread count. The Default is Rayon(0), which is what dot uses internally.

For small matrices, or when you are already inside a parallel region and want to avoid nested threading, pass Parallelism::Serial. For large multiplies on an otherwise idle machine, Parallelism::Rayon(0) lets the engine spread the work. The threading strategy and how the engine picks a thread count are covered in Parallelism in Practice.

Reusing a workspace

The engine needs scratch space to pack blocks of A and B. By default it borrows that space from a thread-local pool, so gemm and dot allocate nothing of their own per call in the steady state. When you run many multiplies in a tight loop and want full control over that buffer, gemm_with takes a &mut Workspace you own and reuse:

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::{Parallelism, Workspace};
use nalgebra::DMatrix;

let mut ws = Workspace::new();
let a = DMatrix::<f64>::from_element(64, 64, 1.0);
let b = DMatrix::<f64>::from_element(64, 64, 2.0);

for _ in 0..1000 {
    let mut c = DMatrix::<f64>::zeros(64, 64);
    gemmkit_nalgebra::gemm_with(&mut ws, 1.0, &a, &b, 0.0, &mut c, Parallelism::default());
}
}

A single Workspace can back multiplies of different shapes across iterations. It grows to fit the largest one it has seen and keeps that capacity. Workspace::new() starts empty and costs nothing until the first call fills it. The _with form exists on all the accumulating entries, including the feature-gated ones. The same reuse pattern carries over to integer, complex, and fused calls.

nalgebra Adapter Advanced Usage

Beyond the real-scalar gemm/gemm_with/dot covered in Using gemmkit with nalgebra, the adapter exposes the engine’s full surface behind Cargo features:

• integer GEMM
• complex GEMM
• fused epilogues
• requantized output
• per-element maps
• prepacked operands
• batching

Every entry keeps the adapter’s core promise. It reads nalgebra’s pointer and strides straight through with no copy. It also mirrors a gemmkit core function of the same name. Cargo gates each family, so you pay only for what you turn on.

Every helper type these entries name comes from the adapter crate itself, so you do not need to name gemmkit to use them. These are Bias, Activation, RequantScale, Requantize, PackedLhs, PackedRhs, Parallelism, Workspace, and Complex (with its c32/c64 aliases).

Integer GEMM

Under int8, gemm_i8 multiplies 2 i8 matrices into an i32 output. The inputs are i8. alpha, beta, and C are i32, because an i8*i8 product needs the wider accumulator. Arithmetic wraps on overflow, the conventional integer-GEMM convention. This is a separate entry from gemm, because the input and output element types differ.

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::{Parallelism, dot_i8, gemm_i8};
use nalgebra::DMatrix;

let a = DMatrix::from_row_slice(2, 3, &[1_i8, 2, 3, 4, 5, 6]);
let b = DMatrix::from_row_slice(3, 2, &[1_i8, 0, 0, 1, 1, 1]);

// dot_i8: A*B into a fresh DMatrix<i32>
let c = dot_i8(&a, &b);

// gemm_i8: scale and accumulate into an i32 output
let mut acc = DMatrix::<i32>::zeros(2, 2);
gemm_i8(1, &a, &b, 0, &mut acc, Parallelism::Serial);
assert_eq!(acc, c);
}

dot_i8(a, b) -> DMatrix<i32> is the allocating convenience form, and gemm_i8_with reuses a caller-owned Workspace for the fixed-cost quantized-inference loop.

Requantized output

Requantization folds the dequantize-scale-round-clamp step of quantized inference into the GEMM, so the m*n i32 accumulator never has to be materialized in full. gemm_i8_requant takes i8 inputs and writes an i8 output. gemm_i8_requant_u8 writes an unsigned u8 output instead (the ONNX QLinearMatMul convention). Both need int8 + epilogue. There is no alpha (it folds into the scale) and no beta (accumulating into a quantized C is ill-defined). The parameters ride in a re-exported Requantize:

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::{Parallelism, RequantScale, Requantize, gemm_i8_requant};
use nalgebra::DMatrix;

let a = DMatrix::from_row_slice(2, 3, &[10_i8, -4, 7, 3, 8, -2]);
let b = DMatrix::from_row_slice(3, 2, &[2_i8, 1, -1, 5, 4, 0]);
let bias = [100_i32, -50]; // per-row, length A.rows

let req = Requantize {
    scale: RequantScale::PerTensor(0.05),
    zero_point: -7,        // in [-128, 127] for the i8 output
    bias: Some(&bias),
};
let mut c = DMatrix::from_element(2, 2, 0_i8);
gemm_i8_requant(&a, &b, req, &mut c, Parallelism::Serial);
}

RequantScale::PerTensor(s) applies 1 scale to every element. RequantScale::PerRow(&[f32]) gives 1 scale per output row (per output channel, the standard per-channel convention), with length A.rows. Every scale must be finite and > 0. The engine adds zero_point as an integer after rounding. It must lie in the output domain: [-128, 127] for gemm_i8_requant, [0, 255] for gemm_i8_requant_u8. The engine adds the optional per-row i32 bias (length A.rows) to the accumulator before scaling.

The adapter validates all of this. It panics with the core engine’s wording on any violation:

• a non-finite or non-positive scale
• a per-row scale or bias of the wrong length
• a zero_point out of range
• a scale or bias slice that overlaps C

Complex GEMM

Under complex, gemm_cplx computes C <- alpha*op(A)*op(B) + beta*C for T = Complex<f32> or Complex<f64>. op(A) is conj(A) when the conj_a flag is set, and op(B) is conj(B) when conj_b is set. The conjugation flags are why complex needs its own entry: they do not fit the homogeneous real-scalar signature. dot_cplx(a, b) is the non-conjugated A*B convenience form. For a conjugated product, use gemm_cplx directly.

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::{Complex, Parallelism, dot_cplx, gemm_cplx};
use nalgebra::DMatrix;

type C = Complex<f64>;
let a = DMatrix::from_element(2, 2, C::new(1.0, 1.0));
let b = DMatrix::from_element(2, 2, C::new(0.0, -1.0));

// plain product
let p = dot_cplx(&a, &b);

// conjugate A, plain B, accumulate into an existing C
let mut acc = DMatrix::from_element(2, 2, C::new(0.0, 0.0));
gemm_cplx(C::new(1.0, 0.0), &a, true, &b, false,
          C::new(0.0, 0.0), &mut acc, Parallelism::Serial);
}

Complex with a fused bias

gemm_cplx_fused (needs complex + epilogue) adds a bias in the same pass as the complex product: C <- alpha*op(A)*op(B) + beta*C + bias. The bias is a re-exported Bias, either Bias::PerRow (length A.rows) or Bias::PerCol (length B.cols). gemmkit adds it verbatim, without conjugation. There is deliberately no activation parameter here: an ordering activation like ReLU is undefined on complex numbers. With bias == None the call is exactly gemm_cplx.

Fused epilogues

The epilogue feature adds gemm_fused, which computes C <- act(alpha*A*B + beta*C + bias) in a single pass over f32/f64. When half is also on, f16/bf16 use the same path. Their epilogue runs in f32, with 1 narrowing store at the end. The optional Bias is PerRow (length A.rows) or PerCol (length B.cols). The optional Activation is Relu or LeakyRelu(slope), applied last. Passing None for both is bit-for-bit identical to plain gemm.

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::{Activation, Bias, Parallelism, gemm_fused};
use nalgebra::DMatrix;

let a = DMatrix::<f32>::from_element(12, 9, 0.5);
let b = DMatrix::<f32>::from_element(9, 7, -0.25);
let bias: Vec<f32> = (0..12).map(|i| 0.5 * i as f32 - 2.0).collect();
let mut c = DMatrix::<f32>::zeros(12, 7);

gemm_fused(1.3, &a, &b, -0.7, &mut c,
           Some(Bias::PerRow(&bias)), Some(Activation::Relu), Parallelism::Serial);
}

The fused pass is not just a convenience. It avoids a second sweep over C, and it avoids the round trip through memory that a separate bias-add and activation would cost. For f32/f64, the result is bit-identical to running gemm and then applying the same bias and activation element by element. You can adopt it without changing numerical results. The design behind it is covered in Fused Epilogues.

Per-element map

gemm_map fits epilogues that a bias-plus-activation shape cannot express. It applies an arbitrary closure to each finished output element: C[r, c] <- f(alpha*A*B + beta*C, r, c). The closure fires exactly once per element, with (r, c) in the user frame of C. T is f32/f64 only. The closure is &(dyn Fn(T, usize, usize) -> T + Sync). It must be Sync to run in parallel, and it can close over data by reference.

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::{Parallelism, gemm_map};
use nalgebra::DMatrix;

let a = DMatrix::<f64>::from_element(8, 6, 0.3);
let b = DMatrix::<f64>::from_element(6, 5, 0.4);
let mut c = DMatrix::<f64>::zeros(8, 5);

// sigmoid, ignoring position
let sigmoid = |v: f64, _r: usize, _c: usize| 1.0 / (1.0 + (-v).exp());
gemm_map(1.0, &a, &b, 0.0, &mut c, &sigmoid, Parallelism::Serial);
}

Prefer gemm_fused for a plain bias or ReLU, because it vectorizes. gemm_map is the general extension point for GELU, sigmoid, clamps, and position-dependent transforms, at the cost of 1 indirect call per output element. Like the fused entry, its f32/f64 result is bit-identical to gemm followed by the same map. gemm_map_with reuses a Workspace.

Prepacked operands

For example, a weight matrix served against a stream of activations has 1 fixed operand across many multiplies. Prepacking it once removes the per-call repack. prepack_rhs(b) -> PackedRhs<T> packs a right operand. gemm_packed_b then consumes the handle in place of B. prepack_lhs/gemm_packed_a do the mirror for a fixed left operand.

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::{Parallelism, gemm_packed_b, prepack_rhs};
use nalgebra::DMatrix;

let weights = DMatrix::<f32>::from_fn(64, 32, |i, j| 0.01 * (i as f32 - j as f32));
let packed = prepack_rhs(&weights); // pack the fixed B once

for step in 0..100 {
    let x = DMatrix::<f32>::from_element(16, 64, step as f32); // an activation batch
    let mut y = DMatrix::<f32>::zeros(16, 32);                 // column-major output
    gemm_packed_b(1.0, &x, &packed, 0.0, &mut y, Parallelism::default());
}
}

There is an orientation constraint. gemm_packed_b needs a column-major-ish C (|col stride| >= |row stride|). A row-major C would force the engine to swap A and B internally. That would invalidate a prepacked RHS, so gemmkit rejects a row-major C. gemm_packed_a is the opposite: it needs a row-major-ish C, and rejects a column-major one. For a C in the wrong orientation, fall back to plain gemm.

Each packed entry has a _with twin for workspace reuse. Under epilogue, the fused twins gemm_packed_b_fused and gemm_packed_a_fused add a bias and activation off the same handle. PackedRhs and PackedLhs expose .rows() and .cols() if you need to re-check dimensions. See Prepacked Operands for the underlying reuse model.

Batched GEMM

nalgebra has no rank-3 array type, so batched GEMM does not take a 3-D tensor the way the ndarray adapter does. Instead, gemm_batched takes the batch as a slice of per-element (&A, &B) input pairs, matched positionally with a slice of &mut C outputs. It runs C_e <- alpha*A_e*B_e + beta*C_e for every element in 1 call, over gemmkit’s pointer-array engine. The alpha, beta, and par arguments are shared by the whole batch.

#![allow(unused)]
fn main() {
use gemmkit_nalgebra::{Parallelism, gemm_batched};
use nalgebra::DMatrix;

let a = DMatrix::from_row_slice(2, 2, &[1.0_f32, 2.0, 3.0, 4.0]);
let b = DMatrix::from_row_slice(2, 2, &[5.0_f32, 6.0, 7.0, 8.0]);
let mut c = vec![DMatrix::<f32>::zeros(2, 2), DMatrix::<f32>::zeros(2, 2)];

let ab = [(&a, &b), (&a, &b)];
gemm_batched(1.0, &ab, 0.0, &mut c, Parallelism::Serial);
assert_eq!(c[0], DMatrix::from_row_slice(2, 2, &[19.0, 22.0, 43.0, 50.0]));
}

Element shapes may differ, so a heterogeneous batch is fine, as long as each element’s own dimensions agree: A_e.cols == B_e.rows, and so on. The shared storage type carries the varying runtime dimensions. DMatrix or dynamic-stride views therefore cover heterogeneous shapes and mixed layouts under 1 type. The input count and the output count must match: ab.len() == c.len(). A mismatch there, or a dimension mismatch in any element, panics.

The batch parallelizes across elements, not within a single element. The dispatcher assigns whole per-element GEMMs to workers, and each worker runs its GEMM serially and cache-hot. As a result, the batch reproduces a plain loop of gemm calls. It also stays reproducible across thread counts. The serial and parallel schedules are bit-identical, because each element always runs wholly on 1 worker.

The C slice is a single storage type. There is no pointer-array analogue of a single shared fused epilogue for it, so unlike the ndarray adapter, there is no gemm_batched_fused here. The batching model is discussed further in Batched GEMM.

Where this sits next to nalgebra’s own multiply

nalgebra already multiplies matrices: &a * &b, a.mul_to(&b, &mut c), and the rest of its operator surface. These return properly typed matrices and integrate with its const-generic dimensions. For an ordinary f32/f64 product, especially of small static matrices, those are the idiomatic choice and there is no reason to reach for this adapter.

The adapter earns its place when you want something nalgebra’s operators do not offer:

• the engine’s runtime SIMD dispatch, which selects the best available instruction set on the machine at run time instead of at compile time
• the i8 -> i32 and requantizing integer paths
• fused bias/activation and per-element epilogues in a single pass
• batched multiplication of many small problems
• prepacked operands for a fixed weight reused across calls

When your matrices already live in nalgebra and you need any of those, reach for the adapter. It gives you the engine’s throughput and features without leaving nalgebra’s types, and without a copy.

Using gemmkit with faer

gemmkit-faer is a thin, zero-copy bridge from faer’s view types to the gemmkit GEMM engine. It accepts a MatRef<'_, T> for each input and a MatMut<'_, T> for the output. It reads the data pointer and the element-unit row and column strides straight out of the view. Then it hands them to gemmkit’s raw engine. The adapter does not transpose, copy, or repack anything on the way in.

faer already stores its strides the way gemmkit’s engine wants them. Because of this, a faer Mat, a transposed view, an offset sub-matrix, and a reversed (negative-stride) view all reach the kernel untouched.

The crate targets faer 0.24. It requires Rust 1.89.

Installation and features

gemmkit-faer re-exports everything its own signatures name. Because of this, a direct gemmkit dependency is not part of the normal setup.

[dependencies]
gemmkit-faer = "0.1"
faer = "0.24"

gemmkit_faer re-exports:

• The Parallelism selector, and the Workspace type every _with variant takes.
• The fused selectors Bias and Activation.
• The prepacked handles PackedLhs and PackedRhs.
• The requantization parameters Requantize and RequantScale.
• The element-type bounds GemmScalar, FusedScalar, MapScalar, and ComplexScalar. Name these when you write a wrapper generic over an entry.
• The element types f16, bf16, Complex, c32, and c64, each under its own feature. This keeps half and num-complex out of your manifest too.
• The tuning module.

Reach tuning through the adapter, not through a separate gemmkit dependency of your own. The knobs are process-global atomics. A second, separately resolved gemmkit gives you a different set of atomics, one the adapter never reads.

Every Cargo feature forwards to the same-named gemmkit feature. So when you enable an element family or a fused entry here, the core turns it on too.

• parallel (default): rayon-based parallelism.
• wasm_threads: threading on wasm32-wasip1-threads. Also enables parallel.
• half: the f16 and bf16 element types, accumulated in f32.
• complex: the c32 and c64 element types.
• int8: i8 inputs into an i32 output.
• epilogue: the fused bias/activation, requantization, and per-element map entries.

The advanced usage page covers the feature-gated families and the fused entries. This page stays on the always-available f32/f64 (plus f16/bf16 under half) surface.

What zero-copy means here

Every entry routes through one small helper. This helper pulls the raw parts out of a MatRef. faer reports strides in element units as isize, negative for a reversed view. This is exactly the shape gemmkit’s unchecked engine expects, so no conversion step exists at all.

#![allow(unused)]
fn main() {
// gemmkit-faer/src/common.rs
pub(crate) fn ref_parts<T>(a: MatRef<'_, T>) -> (usize, usize, isize, isize, *const T) {
    (a.nrows(), a.ncols(), a.row_stride(), a.col_stride(), a.as_ptr())
}
}

The adapter validates the 3 shared dimensions itself. It then forwards the pointers and strides to gemmkit’s _unchecked engine inside a single unsafe block. The safety argument is short. faer’s view types guarantee that the pointer plus strides describe a valid in-bounds layout. The output is a MatMut, an exclusive borrow, so C cannot alias A or B. That is the entire adapter for the plain path.

All of gemmkit’s cache blocking, ISA dispatch, packing, and parallel scheduling live in the core. The core documents them there too. See the architecture chapter for the internals.

gemm and dot

The 2 workhorses are dot, which returns a fresh product, and gemm, which updates an output in place. Both are generic over GemmScalar: f32 and f64 always, plus f16 and bf16 when the half feature is on.

#![allow(unused)]
fn main() {
use faer::Mat;

let a = Mat::from_fn(2, 2, |i, j| [[1.0_f64, 2.0], [3.0, 4.0]][i][j]);
let b = Mat::from_fn(2, 2, |i, j| [[5.0_f64, 6.0], [7.0, 8.0]][i][j]);
// A*B into a fresh column-major Mat
let c = gemmkit_faer::dot(a.as_dyn_stride(), b.as_dyn_stride());
assert_eq!(c[(0, 0)], 19.0);
assert_eq!(c[(1, 1)], 50.0);
}

dot(a, b) computes A*B into a newly allocated column-major Mat. It runs with the default parallelism, Parallelism::Rayon(0), which auto-detects the thread count. Use dot as the one-shot convenience. When you own the output buffer, or want the general update, use gemm instead.

#![allow(unused)]
fn main() {
use faer::Mat;
use gemmkit_faer::{Parallelism, gemm};

let a = Mat::<f64>::from_fn(4, 3, |i, j| (i + j) as f64);
let b = Mat::<f64>::from_fn(3, 5, |i, j| (i as f64) * (j as f64));
let mut c = Mat::<f64>::zeros(4, 5);
// c <- 1.5 * a * b + 2.0 * c, single-threaded
gemm(1.5, a.as_dyn_stride(), b.as_dyn_stride(), 2.0, c.as_dyn_stride_mut(), Parallelism::Serial);
}

gemm(alpha, a, b, beta, c, par) computes C <- alpha*A*B + beta*C in place. With beta == 0, gemm overwrites the prior contents of C and never reads them. This is exactly what dot does internally. With a nonzero beta, the call accumulates onto what C already holds.

The signatures are the ones you see above. The inputs are MatRef<'_, T>, the output is MatMut<'_, T>, and par is a Parallelism. The .as_dyn_stride() and .as_dyn_stride_mut() conversions turn faer’s statically typed strides into the dynamic-stride views the adapter accepts. They cost nothing at runtime.

Layouts that pass through untouched

The adapter only ever reads a pointer and 2 strides. Because of this, any faer view works without a copy or a fallback. A transposed operand is the common row-major-A case. Transposing a column-major matrix yields a view whose row stride is non-unit, and that view goes straight to the kernel.

#![allow(unused)]
fn main() {
// `at` is k x m column-major; `.transpose()` gives an m x k view with a non-unit
// row stride - read straight through, no copy
let a = at.as_dyn_stride().transpose();
let c = gemmkit_faer::dot(a, b.as_dyn_stride());
}

The same holds for an offset sub-matrix. submatrix(...) moves the base pointer and keeps a non-contiguous column stride. It also holds for a reversed view. reverse_rows() and reverse_cols() carry a negative stride.

gemmkit’s unchecked path handles negative strides directly. So a reversed input accumulates correctly under beta, just like any other input. See Matrix Views and Layouts for how the engine treats general strides.

Choosing parallelism

Every entry takes a Parallelism. Parallelism::Serial runs single-threaded. Parallelism::Rayon(n) uses rayon with at most n threads. Rayon(0) auto-detects the thread count.

gemmkit ramps the thread count with the workload, instead of jumping straight to every core. For a fixed machine and a fixed configuration, a call gives reproducible results. Serial and parallel runs also agree bit-for-bit today. That agreement is not a hard guarantee, because the reproducibility contract covers a fixed configuration only, and the worker count is part of that configuration. The Parallelism in Practice guide covers the scheduling model.

Reusing a workspace across calls

gemm allocates its scratch space from a thread-local pool. Every entry also has a _with twin. If you drive many GEMMs in a loop and want to own the scratch buffer explicitly, use the _with twin instead. It takes a &mut Workspace as its first argument and reuses that workspace across calls.

#![allow(unused)]
fn main() {
use gemmkit_faer::{Parallelism, Workspace, gemm_with};

let mut ws = Workspace::new();
for (a, b, mut c) in problems {
    // same result as `gemm`, but the scratch buffer is reused
    gemm_with(&mut ws, 1.0, a, b, 0.0, c.as_dyn_stride_mut(), Parallelism::Rayon(0));
}
}

A single Workspace grows to fit the largest problem it has seen. After that, gemmkit reuses it as-is. This matters most for a stream of similar small-to-medium GEMMs, where allocation would otherwise show up in the profile.

Panic behavior

The adapter checks the 3 shared dimensions before dispatching, and panics on a mismatch. A.cols must equal B.rows. A.rows must equal C.rows. B.cols must equal C.cols. The adapter prefixes each message with gemmkit-faer: and names the 2 conflicting extents, for example gemmkit-faer: A.cols (4) != B.rows (5). These are the only panics on the plain gemm/dot path.

The feature-gated entries add a few more checks: bias length and overlap, requantize parameters, and prepacked-C orientation. Those checks reproduce gemmkit’s own checked-entry wording. Each entry on the advanced usage page lists its own panics.

faer Adapter Advanced Usage

Beyond gemm and dot, the faer adapter mirrors the rest of gemmkit’s surface. This includes the extra element families, the fused epilogues, batched GEMM over a slice, and prepacked operands. Each of these is feature-gated. Each also reads raw pointers and strides straight out of faer’s views. So transposed, sub-matrix, and reversed operands keep working exactly as they do on the plain path. This page walks the families one at a time. It closes with a note on when the adapter earns its place next to faer’s own matmul.

The introductory page covers installation, the zero-copy mechanism, gemm/gemm_with/dot, parallelism, and the workspace pattern. Everything here builds on that page. As on the plain path, every entry also has a _with twin that reuses a caller-owned Workspace.

Integer GEMM (int8)

With the int8 feature, gemm_i8 and dot_i8 take i8 inputs and accumulate into an i32 output. The input and output element types differ. That is why this is a separate entry from gemm, rather than another instance of the generic. faer’s view types are generic over the element. So an i8 MatRef and an i32 MatMut need no special handling.

#![allow(unused)]
fn main() {
use faer::Mat;
use gemmkit_faer::{Parallelism, dot_i8, gemm_i8};

let a = Mat::<i8>::from_fn(16, 12, |i, j| ((i + j) as i8 % 7) - 3);
let b = Mat::<i8>::from_fn(12, 10, |i, j| ((i * 2 + j) as i8 % 5) - 2);
// i8 * i8 accumulated into a fresh Mat<i32>
let c = dot_i8(a.as_dyn_stride(), b.as_dyn_stride());

// Mat::zeros is ComplexField-only, so integer outputs use from_fn
let mut acc = Mat::<i32>::from_fn(16, 10, |_, _| 0);
// c <- 3 * a * b + (-2) * c, all of alpha/beta/C in i32
gemm_i8(3, a.as_dyn_stride(), b.as_dyn_stride(), -2, acc.as_dyn_stride_mut(), Parallelism::Serial);
}

alpha, beta, and C are i32. The arithmetic wraps on overflow. This is the conventional integer-GEMM semantics.

Requantized output (int8 + epilogue)

With both int8 and epilogue, gemm_i8_requant fuses the requantize step into the kernel’s store. i8 inputs multiply into an i32 accumulator. The kernel scales, biases, rounds, and clamps that accumulator down to an i8 output in a single pass. It never materializes the full m*n i32 matrix. gemm_i8_requant_u8 does the same, but clamps to an unsigned u8 output, the ONNX QLinearMatMul-style activation domain.

There is no alpha, because it folds into the scale. There is no beta, because accumulating into a quantized output is ill-defined.

The parameters come in a Requantize. The crate re-exports this type, so you need not depend on gemmkit for it. scale is a RequantScale, either PerTensor(f32) or a PerRow(&[f32]) of per-channel scales. zero_point joins in as an integer after rounding. bias is an optional per-row i32 vector, added to the accumulator before scaling.

#![allow(unused)]
fn main() {
use faer::Mat;
use gemmkit_faer::{Parallelism, RequantScale, Requantize, gemm_i8_requant};

let (m, n) = (17, 13);
let bias: Vec<i32> = (0..m as i32).map(|i| 40 * i - 200).collect();
let mut c = Mat::<i8>::from_fn(m, n, |_, _| 0);
let req = Requantize {
    scale: RequantScale::PerTensor(0.05),
    zero_point: -7,
    bias: Some(&bias),
};
gemm_i8_requant(a.as_dyn_stride(), b.as_dyn_stride(), req, c.as_dyn_stride_mut(), Parallelism::Serial);
}

The output is C[i,j] = clamp(zero_point + round_ne(scale * (sum_k A*B + bias[i])), LO, HI), with round-half-to-even. [LO, HI] is [-128, 127] for the i8 entry and [0, 255] for the u8 entry.

The adapter validates the requantize parameters before dispatch. It reproduces gemmkit’s own checked-entry wording. The checks cover:

• A non-finite or non-positive scale.
• A per-row scale slice of the wrong length, or one that overlaps C.
• A zero_point outside the output domain.
• A bias of the wrong length, or one that overlaps C.

That validation is raw pointer math against C’s byte footprint. The adapter never builds a C slice. This is what lets it forward negative-stride views to the raw engine safely.

Complex GEMM (complex)

With the complex feature, gemm_cplx, gemm_cplx_with, and dot_cplx operate on complex matrices with optional per-operand conjugation. The element type T is Complex<f32> or Complex<f64>.

This is not a separate representation from faer’s own. faer 0.24’s c32 and c64 are type aliases for num_complex::Complex<f32> and num_complex::Complex<f64>. This crate re-exports the same types as Complex, with the same c32/c64 aliases, and constrains its ComplexScalar bound over them. So a faer complex Mat reaches the adapter with no conversion, just like a real one.

gemm_cplx is a separate entry from gemm, because the conjugation flags do not fit the homogeneous surface. It computes C <- alpha*op(A)*op(B) + beta*C, where op(A) = conj(A) when conj_a is set, and op(B) = conj(B) when conj_b is set.

The implementation in cplx.rs pulls the same raw parts as the real path. It threads the 2 bool flags through to gemm_cplx_unchecked. Nothing else differs, so transposed, sub-matrix, and reversed views work identically. dot_cplx is the non-conjugated A*B convenience.

#![allow(unused)]
fn main() {
use faer::Mat;
use gemmkit_faer::{Complex, Parallelism, gemm_cplx};

type C = Complex<f64>;
let a = Mat::<C>::from_fn(12, 9, |i, j| C::new(i as f64, j as f64));
let b = Mat::<C>::from_fn(9, 7, |i, j| C::new((i + j) as f64, 1.0));
let mut c = Mat::<C>::zeros(12, 7);
// C <- alpha * conj(A) * B + beta * C
gemm_cplx(
    C::new(1.3, -0.4),
    a.as_dyn_stride(), true,   // conjugate A
    b.as_dyn_stride(), false,  // leave B
    C::new(0.5, 0.7),
    c.as_dyn_stride_mut(),
    Parallelism::Serial,
);
}

Under complex plus epilogue, there is gemm_cplx_fused. It adds an optional bias in one pass: C <- alpha*op(A)*op(B) + beta*C + bias. The bias is a Bias::PerRow (length A.rows) or a Bias::PerCol (length B.cols). gemmkit adds it verbatim, to every element of that row or column, and never conjugates it.

There is deliberately no activation parameter. An ordering activation such as ReLU is undefined on complex numbers, so the fused complex entry carries a bias only.

Fused bias and activation (epilogue)

With epilogue, gemm_fused computes C <- act(alpha*A*B + beta*C + bias) in a single pass. The optional Bias is PerRow or PerCol. The optional Activation is Relu or LeakyRelu(slope), applied last. Passing None for both gives exactly gemm. The crate re-exports both selectors.

#![allow(unused)]
fn main() {
use gemmkit_faer::{Activation, Bias, Parallelism, gemm_fused};

let bias: Vec<f64> = (0..m).map(|i| 0.5 * i as f64 - 2.0).collect();
// C <- relu(1.3 * A*B - 0.7 * C + rowbias)
gemm_fused(
    1.3, a.as_dyn_stride(), b.as_dyn_stride(), -0.7,
    c.as_dyn_stride_mut(),
    Some(Bias::PerRow(&bias)),
    Some(Activation::Relu),
    Parallelism::Rayon(0),
);
}

For f32/f64, the fused result is bit-identical to plain gemm followed by the same scalar map, for every shape. The epilogue folds into the same kernel’s store without perturbing the accumulation order. Serial and parallel runs also agree bit-for-bit today. That agreement is a property of today’s implementation, not a hard guarantee. The reproducibility contract itself covers only a fixed configuration, and the worker count is part of that configuration.

For f16/bf16 (under half), the fused result is more precise instead of identical. A separate gemm() call followed by a narrow map rounds to the narrow type, widens back, and rounds again. The fused path skips that extra rounding. The bias and slope widen exactly to f32, the epilogue applies in f32, and the result rounds once to the narrow output. So for f16/bf16 the fused result is more precise, though it is not bit-identical to gemm followed by a narrow map. The f32/f64 bitwise guarantee above does not extend to these narrow types. Serial and parallel runs still agree bit-for-bit for these types too, under the same fixed-configuration reproducibility contract. The Fused Epilogues guide has the full contract.

For an arbitrary per-element function, there is gemm_map (f32/f64 only): C[r,c] <- f(alpha*A*B + beta*C, r, c). The closure runs once per output element, at its final value, with (r, c) in the user frame of C.

Use gemm_map for GELU, sigmoid, clamps, or position-dependent transforms. Prefer gemm_fused instead for a plain bias or ReLU, because it vectorizes. gemm_map pays one indirect call per element.

Batched GEMM

faer has no rank-3 array type, so gemmkit-faer expresses batched GEMM over slices instead. gemm_batched takes a &[(MatRef, MatRef)] of per-element (A, B) inputs, paired positionally with a &mut [MatMut] of C outputs. All elements share one alpha, beta, and Parallelism.

gemmkit’s pointer-array engine parallelizes the batch across elements. Its scheduler assigns whole GEMMs to workers. Each worker runs its GEMM serially and stays cache-hot for it.

#![allow(unused)]
fn main() {
use faer::Mat;
use gemmkit_faer::{Parallelism, gemm_batched};

let a = Mat::from_fn(2, 2, |i, j| [[1.0_f64, 2.0], [3.0, 4.0]][i][j]);
let b = Mat::from_fn(2, 2, |i, j| [[5.0_f64, 6.0], [7.0, 8.0]][i][j]);
let mut c0 = Mat::<f64>::zeros(2, 2);
let mut c1 = Mat::<f64>::zeros(2, 2);
let ab = [
    (a.as_dyn_stride(), b.as_dyn_stride()),
    (a.as_dyn_stride(), b.as_dyn_stride()),
];
let mut c = [c0.as_dyn_stride_mut(), c1.as_dyn_stride_mut()];
gemm_batched(1.0, &ab, 0.0, &mut c, Parallelism::Serial);
}

Element shapes may differ, a heterogeneous batch, as long as each element’s own dimensions agree. The call panics if the input and output counts disagree. It also panics if any element’s dimensions are inconsistent, naming the offending element index.

Each element re-dispatches through the full engine. So the batch reproduces a plain loop of gemm calls. It is deterministic across thread counts, because each element runs wholly on one worker. For the same reason, serial and batch-parallel output are bit-identical.

There is no batched fused entry here. The ndarray adapter offers a shared-epilogue batched form, but it has no pointer-array analogue in the core. See Batched GEMM for the scheduling policy.

Prepacked operands

When one operand stays fixed across many calls, for example weights against a stream of activations, pre-pack it once and skip the per-call repack. prepack_rhs turns a B into a reusable PackedRhs, consumed by gemm_packed_b. prepack_lhs turns an A into a PackedLhs, consumed by gemm_packed_a. The crate re-exports both handles.

#![allow(unused)]
fn main() {
use gemmkit_faer::{Parallelism, gemm_packed_b, prepack_rhs};

let packed = prepack_rhs(weights.as_dyn_stride()); // pack the fixed B once
for (act, mut out) in stream {
    // out must be column-major-ish (|col stride| >= |row stride|)
    gemm_packed_b(1.0, act.as_dyn_stride(), &packed, 0.0, out.as_dyn_stride_mut(), Parallelism::Rayon(0));
}
}

The one constraint is output orientation. A prepacked B fixes the operand roles, so gemm_packed_b needs a column-major-ish C (|col stride| >= |row stride|). A row-major C would swap the A/B roles and invalidate the packed RHS, so gemmkit rejects it. Symmetrically, gemm_packed_a needs a row-major-ish C. For a mismatched output layout, fall back to plain gemm.

Under epilogue, the prepacked entries have fused twins: gemm_packed_b_fused and gemm_packed_a_fused. Each takes the same Bias/Activation as gemm_fused, off the same handle. The Prepacked Operands guide explains the reuse model.

When to reach for this adapter

faer ships its own matmul. For a plain f32/f64 product of 2 faer matrices, use that instead. This adapter earns its place when you need something the core faer operator does not offer, on faer’s own types, without leaving the faer ecosystem:

• Extra element families: i8 -> i32 integer GEMM, and requantization that fuses straight down to an i8 or u8 output.
• Fused epilogues: the kernel computes bias and activation, or an arbitrary per-element closure, in the same pass as the product, not as a second sweep over C.
• Prepacking across calls: pack a fixed weight matrix once, then reuse it over a long inference loop.
• A shared tuning surface: all 3 gemmkit adapters sit on the same engine, so one GEMMKIT_* environment profile from gemmkit-tune applies to all of them. See Tuning Knobs for the knob surface.

If none of those apply, use faer’s built-in matmul instead. It is the simpler choice. This adapter supplements it. It does not replace it.

Tuning with gemmkit-tune

gemmkit-tune is a small command-line autotuner. You run it once on the machine that will run your gemmkit workload. It sweeps gemmkit’s runtime GEMMKIT_* knobs. For each knob, it measures a representative set of matrix shapes, then writes a shell profile of export GEMMKIT_*=... lines. Source that file before you launch a gemmkit binary. The already-built binary is then retuned for the host, with no recompile, no code change, and no dependency on how the binary was distributed.

Why it exists

gemmkit’s compiled-in defaults are not arbitrary. Every threshold in gemmkit::tuning was hand-calibrated against real measurements on one reference machine. The value that is optimal on that machine encodes its cache sizes, its core count, and its ratio of DRAM bandwidth to compute.

A different CPU has a different L2, a different core count, and a different memory bandwidth. Its crossovers land in different places, for example:

• the k at which packing starts to pay
• the problem size at which a shared pre-pass beats per-worker packing
• the byte floor below which a bandwidth-bound gemv should stay single-threaded

gemmkit-tune re-discovers those crossovers on the silicon in front of it and pins them.

There is a useful corollary. Run the tool on the reference machine and it re-selects essentially the shipped defaults. The report says every knob kept its default. That is the correct outcome, not a disappointment. It validates the tool.

The payoff comes on a machine unlike the reference machine. Examples include a laptop, a cloud instance sharing a socket, a wide Graviton, and an Apple part with a shared cluster-L2 and no L3. The further the deploy host is from the reference machine, the more there is to win.

The sweep contains no randomness. Given the same machine and the same flags, it produces the same profile every time. A profile is therefore a reproducible artifact you can commit alongside a deployment.

Why the deploy host, and never build.rs

The knobs are calibrated against the CPU that executes the work, so the tool has to run there. 2 consequences follow.

First, do not run it in a build.rs. The build host is usually not the deploy host. You compile on a CI runner or a developer laptop, then ship the binary to something else entirely. A build.rs autotune would measure the builder instead, then bake the builder’s crossovers into a binary that runs somewhere with a different cache hierarchy. A cross-compiled build cannot even execute the target’s code. The whole point is to measure the real silicon, so the tool must run on it.

Second, the knobs are plain runtime env vars, read once per process at startup. You do not need to rebuild to apply a profile. The same shipped binary reads whatever GEMMKIT_* values are in its environment. Tune the host, source the profile, and launch. The engine reconfigures itself.

Install and a first run

Install the binary and run it on the target machine:

cargo install gemmkit-tune
gemmkit-tune

A full run takes a minute or two and prints a report as it goes. By default it writes gemmkit-tune.env in the current directory. Source that file in the shell that launches your application:

source gemmkit-tune.env
./your-gemmkit-app

That is the whole workflow. Everything below is refinement. It covers bounding the run, matching it to how you deploy, and knowing when to do it again.

The flags

gemmkit-tune takes no positional arguments. All behavior is on 5 flags.

--threads <n>

Tune for this worker count. Every parallel probe runs under Parallelism::Rayon(n), so the scheduling knobs (the oversample factors and the auto worker-count ramp) are optimized for exactly that width. The default is the machine’s available parallelism. It is capped at the machine width. You cannot tune for more workers than the box physically has. The stamped worker count is always truthful.

The rule of thumb: run the sweep with the same worker count your application will actually use. If your app pins Parallelism::Rayon(8), pass --threads 8. A profile tuned for 32 workers can pick a different scheduling grain than one tuned for 8, and the mismatch costs you.

--time-budget <dur>

Cap the sweep and coarsen it to fit. It accepts 30s, 2m, 1h, or a bare number of seconds. Under a budget the tool takes fewer timing repetitions per estimate: 7 by default, 5 under 90 seconds, and 3 under 30 seconds. Once the deadline passes, it stops sweeping and lists the remaining knobs as skipped for “time budget exhausted”.

Use it when install time must be bounded. Leave it off for the most reliable profile. If the budget is so small that not even one knob gets measured, the report says so and tells you to raise it.

--large-matrices <GiB>

Opt into the 2 memory-heavy probes, GEMMKIT_K_STREAM_MAX and GEMMKIT_SHARED_LHS_MNK, with the given GiB figure as the budget for the giant gemv matrices. These knobs only bite in an expensive regime. One needs a gemv output that spills the last-level cache, which means multi-gigabyte matrices. The other needs a very high-FLOP shape above the shared-pre-pass crossover. Both are off by default.

If the budget you pass cannot hold the required probe, the tool skips that knob cleanly and prints the exact GiB figure to re-run with. The GEMMKIT_K_STREAM_MAX probe is 64-bit only. GEMMKIT_SHARED_LHS_MNK still sweeps on 32-bit. Start with 4 or 8 and follow the advice if it asks for more. See Inside the Sweep for what these 2 probes actually do.

--out <path>

Write the profile somewhere other than ./gemmkit-tune.env.

--dry-run

Run the full sweep and print the report, but write no profile. This is good for previewing what a machine would choose before you commit a file. -h / --help prints usage.

Anatomy of the emitted profile

The file is a header comment, followed by one export line per swept knob, followed by a footer listing what was not swept. It looks like this. The values depend entirely on the host:

# gemmkit-tune profile. Source this before you run a gemmkit app: `source <this file>`
# generated 2026-07-19 14:12:03 UTC by gemmkit-tune 0.1.2
# host: 16 logical cores; L1d 32 KiB, L2 1024 KiB, L3 32 MiB; page 4 KiB
# tuned for 16 worker(s)

export GEMMKIT_MC_REG_PANELS=8  # default (1.00x)
export GEMMKIT_LHS_PACK_THRESHOLD=256  # tuned (1.07x)
export GEMMKIT_PAR_MNK_PER_WORKER=4000000  # tuned (1.03x)

# not swept on this host:
#   GEMMKIT_PARALLEL_THRESHOLD: serial/parallel break-even is strongly shape-dependent ...
#   GEMMKIT_DEEP_KC_BYTES: narrow-only (f16/bf16 deep-contraction twin); no narrow probe here ...

The header records the host it was tuned on. It lists the logical core count, the 3 cache levels, the page size, the worker count, the tool version, and a UTC timestamp. That stamp tells you, months later, whether a profile still matches the box you are looking at.

Each export line carries a trailing comment. The comment marks the line default when the winner equals the shipped default, or tuned when the winner moved. It also names the measured speedup over the default. A knob that kept its default is still written. The profile is therefore a complete, self-documenting record of the decision, not just the deltas.

The values are always raw integers. An “unbounded” winner is written as its numeric value, never as a MAX alias, because gemmkit’s env parser reads a plain decimal integer. A malformed GEMMKIT_* value is not fatal. gemmkit warns once on stderr and falls back to the compiled default, so a hand-edited typo degrades to the default instead of crashing.

Deploying the profile

gemmkit reads each GEMMKIT_* variable once, on first access, then caches it for the life of the process. The profile must therefore be in the environment before the first GEMM call. Sourcing it before launch guarantees exactly that. There are 3 common ways to arrange it.

Shell profile or launch script. This is the direct case. Run source gemmkit-tune.env in the same shell that runs the binary, or in the service’s launch script. Full shell semantics apply, so the file drops in unchanged.

Container entrypoint. Bake the profile into the image and source it in the entrypoint before you exec your app. Every container then starts pre-tuned. Tune on a host that matches the container’s runtime hardware, not the image builder.

systemd EnvironmentFile. This works, with one caveat. systemd’s EnvironmentFile parser wants bare NAME=value lines. It does not understand the export keyword or the trailing # tuned (...) comments. Convert the profile first, for example with grep '^export' gemmkit-tune.env | sed -e 's/^export //' -e 's/[[:space:]]*#.*$//' > gemmkit.env, then point EnvironmentFile= at the result. The #-comment header is fine to leave in. Only the assignment lines need this transform.

One note on precedence. A GEMMKIT_* env var is overridden by a programmatic tuning::set_* call in the application. If your app tunes a knob in code, the profile does not change that knob. This is deliberate: self-tuning code wins over a deployment profile. An app that wants the profile to apply simply does not call the setters. See Tuning Knobs for the full precedence order.

Run in a clean environment

Any GEMMKIT_* variable already set in the tuning shell skews the sweep, because gemmkit reads it while measuring the baseline. The tool neutralizes the knobs it sweeps and warns you about any GEMMKIT_* variable it finds set. The reliable move is still to tune from a shell with none of them present. Do not source a previous gemmkit-tune.env and then re-run the tool in the same shell. That is exactly the polluted baseline the warning is about.

When to retune

Retune when the thing the profile was stamped for changes. That means a different deploy machine, since a different CPU brings different cache sizes, which is the whole reason the tool exists. It also means a different worker count: a profile tuned for 8 workers is not the right one for 32.

A gemmkit or gemmkit-tune version bump can also add knobs, so regenerate the profile after upgrading. A profile that no longer matches its header stamp is a profile to throw away and regenerate.

To understand what the sweep actually measures, how it scores candidates, and how to sanity-check a profile against your own workload, read Inside the Sweep.

Inside the Sweep

Tuning with gemmkit-tune covers how to run the tool. This page is the mechanism behind the sweep. It covers:

• what the sweep measures
• how it decides a winner
• why it is biased toward the shipped default
• how to check that a profile actually helps you

One knob at a time

The sweep is a set of independent one-dimensional searches, not a joint optimization. For each knob, every other knob stays at its default. The tool measures the knob’s candidate values back to back and chooses a winner. Then it restores the knob to its default before it sweeps the next one. Each knob is therefore evaluated against an otherwise-default engine.

This is a deliberate simplification. A full joint search over about 30 knobs is combinatorially hopeless, and it would be dominated by noise. The crossovers these knobs gate are, by design, individually meaningful, so a one-dimensional search suits them.

The cost is that the sweep does not explore cross-knob interactions. This is an acceptable trade, because the defaults already sit at a good joint operating point. The tool’s job is only to move individual crossovers to where the host puts them.

The tool measures candidates in a fixed order with no randomness. The default value goes first, since it is the tie-break incumbent, then each distinct extra candidate follows. It rebuilds the buffers for every shape with the same seeds. An A/B comparison between 2 candidate values therefore sees byte-identical inputs, and any machine drift cancels out.

The sweep table stays in lockstep with the engine

gemmkit enumerates its knobs in one place: gemmkit::tuning::knob_env_names(). This machine-readable registry is the single source of truth for every GEMMKIT_* name. The tuner classifies each knob as either TUNED, meaning it has a real sweep, or NEVER_TUNED, with a reason. A test asserts that these 2 lists partition knob_env_names() exactly. No knob is missing, and no entry is stale.

The practical guarantee is direct: a knob added to gemmkit cannot silently escape the autotuner. The build fails until someone writes a sweep for it, or records why it is deliberately left alone. So when you read the tool’s list of swept knobs, you are reading a list the compiler keeps honest against the engine.

What is measured, and in what unit

Each candidate’s score is a throughput. A GEMM, i8, or batched probe is scored in GFLOP/s, computed as 2*m*k*n per call, times the batch count for a batched probe. A gemv probe is scored in GB/s instead, because a matrix-times-vector is bandwidth-bound, and the bytes moved is the honest figure of merit there.

A single-shape estimate is deliberately robust. First, the tool warms up the probe closure a few times. Then it auto-sizes an iteration count so one timed batch runs for about 50 ms. It times several such batches and reports the median rate, along with the observed min and max. The min and max are not cosmetic. They record the run-to-run spread, and the winner logic uses that spread to stay honest under noise.

Scoring: geometric mean over a probe-shape set

A knob is never judged on one shape. Each knob carries a small set of probe shapes. The tool chooses these shapes so the knob actually binds and so its crossover is bracketed on both sides. A candidate’s score is then the geometric mean of its per-shape median throughputs.

The geometric mean gives every shape equal weight, regardless of its absolute size. One big shape therefore cannot flatter a value that only helps that shape. A winner has to be a broad improvement across the whole set. The worst shape’s spread carries through the geomean, so the noise gate stays conservative across the whole set instead of trusting only the calmest shape.

The probes are picked per knob to make the knob bind. A few examples:

The tie-break is default-biased and noise-aware

Picking the highest geomean would be wrong. On a noisy machine, a 1% edge is usually luck. The winner logic instead starts at the default. It upgrades to a candidate only when that candidate’s geomean beats the current best by more than the larger of the 2 candidates’ measured spreads. Run-to-run noise cannot clear that bar by construction, so it can never rewrite a knob. An exact tie keeps the default.

There is a further margin for the “auto” knobs, whose default is 0. These knobs derive their value from the machine, such as LLC size, core count, or page size. A fixed candidate must beat auto by an extra 5% beyond noise. These auto derivations adapt to shapes the probe set does not cover. A fixed number that wins by a hair on the probes is not worth trading that adaptivity for.

This default bias is the right call under noise. The default is a known-good, deliberately chosen value, and the tool often runs unattended on a machine nobody is watching. The asymmetric bar means the worst case is that the tool just reproduces the defaults. It never regresses you into a measurement artifact.

The sweep also has no RNG anywhere, so a run is safe to trust. At worst it does nothing. When it does move a knob, a real and repeatable improvement cleared the noise.

How the time budget caps and coarsens the sweep

--time-budget acts in 2 ways.

First, it coarsens each estimate up front: 7 timing repetitions with no budget, 5 under 90 seconds, and 3 under 30 seconds. This trades a little measurement stability for speed.

Second, it enforces a hard deadline. Before each knob, the tool checks the clock. Once the deadline has passed, it stops starting new sweeps and records every remaining knob as skipped for “time budget exhausted”.

A tight budget therefore both blurs the measurements it does take and drops knobs off the tail. With no budget, the sweep runs to completion at full repetitions.

Which knobs are skipped, and why

Some knobs are never swept. The report and the profile footer say why for each:

• PARALLEL_THRESHOLD: the serial/parallel break-even is strongly shape-dependent. A single m*n*k scalar cannot fit every aspect ratio, so the tool keeps the calibrated cross-shape default instead of auto-fitting it. Contrast GEMV_THRESHOLD, which is a clean binary on/off decision and is swept.
• DEEP_KC_BYTES: this gates the f16/bf16 deep-contraction twin, and the tuner runs no narrow-type probe. Its auto default derives from L2, a machine property. Override it directly if you need to retune the narrow deep-k engage point.
• PREFETCH_MIN_BYTES: this gates the driver’s C-tile prefetch. Its auto default derives from the detected LLC, a machine property, and probing the crossover would need a beyond-LLC working set on every candidate. Override it directly to retune the engage point (usize::MAX disables the prefetch, and 1 forces it on).

Other knobs are inert on the current target and get skipped for that reason. SEQ_INTERNAL_BYTES_PER_WORKER is read only by the aarch64 batched-split planner. It is swept there, and it is inert and skipped on x86. I8_VNNI_MIN_PAR_MNK gates the x86 VNNI small-parallel fallback, which no other target’s i8 kernel has. NC_NO_L3_PANELS is consulted only on a machine with no L3. It is swept there, and it is inert and skipped on an L3 host.

The 2 heavy knobs are skipped unless you pass --large-matrices.

What –large-matrices unlocks

2 knobs only matter in a regime that is expensive to reproduce, so they are opt-in behind a memory budget.

K_STREAM_MAX caps how far the axpy-gemv output stays register-blocked. It only wins once the output is clearly DRAM-bound. Its probe therefore fixes the output at about twice the last-level cache. A 1x-LLC output sits on the cache boundary and measures nothing decisive, so the probe avoids that size. The probe then sweeps k around the calibrated ceiling.

That output size is fixed, not budget-scaled, so reaching it takes multi-gigabyte matrices. If the budget you passed cannot hold the largest probe, the tool skips the knob. It prints the GiB figure you need to re-run with, rounded up. On a 32-bit target, it skips the knob outright, because the matrices do not fit the address space at all.

SHARED_LHS_MNK gates the shared-LHS pre-pass. This pre-pass removes redundant per-worker A-packing, but it adds a fork-join barrier. So it only pays off above a large m*n*k value (about 8e9 on x86). Its probes use tall, high-FLOP shapes above that crossover.

The tool neutralizes both knobs during the ordinary sweeps, whether or not it is sweeping them itself. This way a stale env value cannot skew a baseline that reads them.

Reading the terminal report

The report opens with a one-row-per-knob summary table. Its columns are knob, unit, shape count, default, winner, speedup, and a result column that reads keeps default or -> <value>, with moved knobs highlighted.

Below the table, a candidate-detail block prints the full sweep landscape for each knob. It lists every candidate’s geomean median. A leading mark flags the default value, and a separate mark flags the winner, so you can see how flat or sharp the optimum was.

After that comes the skipped list with reasons, then a footer. The footer counts how many knobs were swept, how many moved off default, and how many were skipped. On the reference machine, the footer notes that all knobs kept their defaults, which is expected, and the profile reproduces them.

Sanity-checking a profile

The sweep measures synthetic, roughly-square probes. That is the right choice for finding a machine’s crossovers. Your workload has its own shapes, though, so confirm the win transfers before you trust a profile in production.

There are 2 ways to check. The direct one: time your own application with and without gemmkit-tune.env sourced, on the deploy host, and compare. The reproducible one: run gemmkit’s criterion benches, which cover 5 headline groups (sgemm, dtypes, gemv, prepacked, batched) under a saved baseline:

cargo bench -p gemmkit -- --save-baseline stock
source gemmkit-tune.env
cargo bench -p gemmkit -- --baseline stock

If a knob moved and something you care about regressed, the profile is a plain text file. Delete or comment out that one export line and keep the rest. The header stamp and per-line tuned/default tags make it easy to see which line to touch.

Design Goals and the Big Picture

gemmkit is a pure-Rust GEMM engine. It computes C <- alpha*A*B + beta*C over &[T] slices with explicit strides, or over raw pointers with isize strides. At runtime it selects the best instruction set the machine offers. On x86-64 that is AVX-512 or FMA/AVX2, with dedicated VNNI and BF16 dot kernels. On aarch64 it is NEON. On wasm32 it is simd128. A portable scalar path runs everywhere else. The workspace uses edition 2024, rust-version 1.89, and the MIT OR Apache-2.0 license.

3 kinds of callers shape the API surface. Application code uses the safe slice entries, such as gemm, gemm_fused, and gemm_i8. These entries validate everything before any unsafe work runs. Linear algebra libraries use the *_unchecked tier. This group includes the shipped ndarray, nalgebra, and faer adapters, plus anything built the same way. The *_unchecked tier trusts the caller’s own invariants and accepts layouts the safe tier cannot express. Constrained deployments get a core that builds #![no_std] with zero mandatory dependencies, down to wasm32 with compile-time SIMD. Everything else in this chapter follows from 4 design tenets. ARCHITECTURE.md states them compactly under “Goals and constraints”. This chapter expands each one with the reasoning behind it.

Safety at the boundary

The checked entries run validate_gemm_views (gemmkit/src/api.rs) before touching any unsafe code. Its panic catalog is deliberately exhaustive:

• Shape mismatch: gemmkit checks A.cols != B.rows, A.rows != C.rows, and B.cols != C.cols. Each panic message names the 2 numbers that disagree.
• A view addressing outside its slice: for A, B, and C, gemmkit computes the highest offset each view’s strides can reach (extent). It checks that offset against the slice length. A view that needs more elements than its slice holds panics with the exact shortfall.
• Negative strides: the safe tier rejects these, with a message that points to gemm_unchecked. A &[T] view with a negative stride would address below the slice start. The safe extent math cannot vouch for that address.
• A self-aliasing output: a stride on C can map 2 distinct (i, j) pairs to the same offset. A zero stride is the common case. This is fine on A or B, since a broadcast input is only read. On C it panics, because the parallel driver assumes output tiles do not overlap. Writing through such a view would create a data race that entirely safe code could reach.
• C overlapping A or B: gemmkit checks this as byte ranges. This stays exact even when C (i32) and A/B (i8) have different element sizes. The fused entries also check bias length (PerRow = m, PerCol = n) and bias and C disjointness.
• A problem too large to size: broadcast strides allow logical dimensions near isize::MAX. The internal pack-buffer sizing can then overflow usize. Every such product panics, fail-closed, at the element-to-byte chokepoint (Workspace::regions), instead of wrapping and under-allocating.

The panic wording is itself a tested contract. The correctness suite asserts the exact strings. Changing an error message is therefore a deliberate, visible act.

The *_unchecked tier exists because this validation is only meaningful at one boundary. The adapters pull pointers and strides straight out of ndarray, nalgebra, and faer types. Those types already guarantee validity through their own invariants, so re-checking would be pure overhead. The slice-based checks could not even express what the adapters need. For example, a reversed ndarray view has a negative stride and a base pointer in the middle of its allocation. Both are legal and sound for the raw engine. Safety is therefore paid exactly once, either by gemmkit’s validator or by the caller’s type system, never both.

The unchecked entries are ordinary unsafe fns with documented contracts. The guide covers them in The Unchecked Tier.

Reproducible, not bitwise, parallel results

gemmkit promises reproducible parallel results. For a fixed input, a fixed environment, and a fixed configuration, the output does not depend on the worker count. Three mechanisms carry this promise.

The first mechanism is the blocking sizes. KC and NC come from the cache model alone and never depend on the thread count. MC only ever changes by an MR-aligned regroup. Every run therefore reduces each output element in the same fixed order.

The second mechanism is the reduction order. One worker reduces each output element start to finish, over the full depth. The engine never splits a reduction across workers.

The third mechanism is demand-driven scheduling. Packed bytes do not depend on who packs them, so any worker can take any tile. Which worker computes a tile varies from run to run. The result never does.

Just as important is what gemmkit does not promise. Bitwise serial-versus-parallel identity is not part of the contract. It happens to hold on the driver paths today, because serial and parallel runs execute the same kernel over the same blocking. Nothing pins this fact in place. Bitwise identity across configurations is explicitly out too. Change a tuning knob, and the blocking, and so the floating-point summation order, may legitimately change.

Bitwise identity across kernels of the same type is out as well. The bf16 vdpbf16ps dot kernel reshapes the accumulation rounding relative to the widen-and-FMA path. gemmkit holds that kernel to a tolerance, not to exact equality.

Why draw the line there? A permanent promise of bitwise serial-versus-parallel identity would forbid useful engineering. It would rule out dot-product instructions that fuse depth pairs. It would also rule out a blocking choice that takes parallelism into account. Such a promise would buy nothing a user can rely on across machines or library versions anyway.

Parallelism-aware blocking is no longer hypothetical. The driver already has a job-depth floor that shrinks MC with the worker count, to keep the parallel job list deep enough. This stays bitwise-reproducible precisely because the weaker contract leaves room for it. MC still stays an MR multiple, so the microtile set and every element’s KC-shaped accumulation order stay unchanged.

Reproducibility under a fixed configuration is the property tests can assert. Deployments can depend on it, and the engine can keep it while it evolves. Where a path can promise more cheaply, it does. gemv partitions output rows and is bit-identical across worker counts. The i8 integer path uses exact arithmetic, so its VNNI dot kernel is bit-identical to the widen kernel.

No macros, no transmute at the variation points

The engine varies along 3 axes: instruction set, element type, and operation family. Each axis is an ordinary trait. Simd and SimdOps cover the ISA. Scalar covers the element type. KernelFamily covers the operation family.

Dispatch slots are typed function pointers, cached in OnceLocks. Microtile geometry is a pair of const generics, chosen at the dispatch site. Here is what a “kernel variant” actually looks like, from gemmkit/src/dispatch/float.rs:

#![allow(unused)]
fn main() {
unsafe fn gemm_f32_fma(t: Task<f32>, par: Parallelism, ws: &mut Workspace) {
    // MR = 2*8 = 16, NR = 6 -> 12 acc + 2 lhs + 1 rhs = 15 of 16 YMM
    unsafe { run_typed::<f32, Fma, 2, 6>(Fma, t, par, ws) }
}

unsafe fn gemm_f32_avx512f(t: Task<f32>, par: Parallelism, ws: &mut Workspace) {
    // MR = 2*16 = 32, NR = 12 -> 24 acc + 2 lhs + 1 rhs = 27 of 32 ZMM
    unsafe { run_typed::<f32, Avx512F, 2, 12>(Avx512F, t, par, ws) }
}
}

That is the entire per-(type, ISA) surface. Each variant is one wrapper that names a token and a tile. The alternative, macro-stamped or hand-copied per-ISA kernels in the C BLAS tradition, was rejected for its cost to review and to extend.

Traits and const generics leave exactly one generic microkernel to read, step through, and fix. A scheduling improvement lands once, and every ISA inherits it. The compiler type-checks every monomorphization. The OnceLock slots hold typed function pointers, not type-erased ones, so a signature drift is a compile error, not a latent transmute bug.

Extension follows the same shape. A new ISA needs a zero-sized token, its SimdOps impls, and one arm per selection ladder. A new element type needs a Scalar impl, a family (or a reuse through the widen/narrow seam), and a dispatch slot. The driver, packing, and blocking never change. A test (gemmkit/tests/open_closed.rs) enforces this by driving the driver with a second, trivial family. 2 follow-up pages walk each seam in detail: SIMD Tokens and ISA Dispatch and Scalars and Kernel Families.

no_std and a zero-mandatory-dependency core

With default features off, the core crate builds #![no_std]. It needs only core and alloc, and depends on nothing else. Every optional feature pulls in at most one crate:

Without std, compile-time target features replace runtime CPU detection. The env knobs turn off, though the programmatic tuning::set_* setters still work, since they are plain atomics. A per-call workspace replaces the thread-local pool.

A kernel this low in the stack should not force a dependency policy on its hosts. An embedded or wasm deployment gets the same driver, the same families, and the same reproducibility contract as a desktop build. It loses only the machinery that genuinely needs an OS. The practical how-to lives in no_std and WebAssembly.

The workspace map

5 crates release in lockstep at version 0.1.2, plus a fuzzing crate that deliberately sits in its own workspace root:

The adapters are thin by design. Each one pulls the matrix pointer and strides straight out of the host library’s native view. That view may be C-order, F-order, general strides, or reversed strides, and the adapter never copies data. Each adapter then forwards to the *_unchecked engine, relying on the host type’s own invariants for the validation the safe tier would otherwise do. Each adapter also forwards the same-named Cargo features (parallel, wasm_threads, half, complex, int8, epilogue) to gemmkit, so the feature story stays identical everywhere. The adapter chapters cover their full surfaces. See ndarray, nalgebra, and faer.

gemmkit-tune is the out-of-process calibrator. Every heuristic threshold in the engine is a runtime knob (see Tuning Knobs). The compiled defaults were calibrated on one machine, so the tuner exists to redo that calibration on yours.

Run the tuner binary once on the deploy host. It sweeps each knob over a set of probe shapes, then writes a gemmkit-tune.env profile of export GEMMKIT_*=... lines. Source that file before you launch your application. This needs no recompile and no build-time coupling. The only contract between the tuner and the library is the documented env-var surface. The tuning::knob_env_names registry keeps that contract honest, since the tuner’s sweep table is checked against it. The gemmkit-tune chapter has the practical guide.

The fuzz crate sits outside the workspace on purpose. cargo-fuzz needs nightly, for build-std and AddressSanitizer. Excluding the fuzz crate keeps cargo test --workspace and the MSRV build on stable.

This chapter and ARCHITECTURE.md

The repository’s ARCHITECTURE.md is the compact map. It gives the layer table, the call path, the seams, and one section per subsystem. It is written for a reader who has the code open in another pane. This book chapter is the guided tour of the same material. It uses the same layer labels and the same file references, but leaves room for the reasoning, the rejected alternatives, and worked examples. When the two disagree, the code wins, and both documents have a bug.

Read The Layer Stack next for the structure. Then read Life of a GEMM Call for the motion.

The Layer Stack

Every module in the core crate opens by declaring its place in a stack. api.rs says “Public core API (layer L8a)”. driver.rs says “The generic GEMM driver (layer L5)”, and so on down to simd.rs at L0. The labels are not decoration. They write the crate’s dependency discipline where a reader cannot miss it. The map below lists the modules in dependency order, which is what makes the downward claim checkable. This page walks the stack from the bottom up. By the time it reaches the public API, every word the API uses has already been defined. The next page, Life of a GEMM Call, traverses the same stack in the other direction, following one call.

L8a  api        safe slice entries, *_with, *_unchecked; MatRef/MatMut
L7   dispatch   runtime ISA selection, one memoized fn pointer per type
L6   special    gemv, small-k, small-m,n, batched reroutes
L5   driver     the generic 5-loop blocked GEMM, one for all families
L4   kernel     KernelFamily seam (float/mixed/int/complex) + Epilogue
L3   cache      topology detection + BLIS analytical blocking
L2   parallel   worker-count resolution, JobCursor work distribution
L1   pack       micropanel packing primitives
L0   simd       ISA tokens + SimdOps vocabulary;  scalar: Scalar/Acc types
     ---        cross-cutting: tuning (GEMMKIT_* knobs), workspace (buffers)

2 placements are worth calling out, because they are what keep every arrow pointing down. parallel sits low, below the kernel families, the cache model, and the driver, because it is self-contained worker vocabulary. The Parallelism policy enum, the Ptr Send-pointer wrapper, and the JobCursor depend only on tuning. kernel, driver, special, and dispatch all reach down to them. pack sits below kernel because the families’ pack hooks build on the packing primitives, not the other way around.

L0: the vocabulary, scalar.rs and simd.rs

The bottom layer defines what the rest of the crate is allowed to talk about. gemmkit/src/scalar.rs holds the data-type seam, and it is deliberately tiny:

#![allow(unused)]
fn main() {
pub trait Scalar: Copy + Send + Sync + PartialEq + 'static {
    /// The type in which products are accumulated. `Self` for `f32`/`f64`
    type Acc: Scalar<Acc = Self::Acc>;
    /// The additive identity
    const ZERO: Self;
    /// The multiplicative identity
    const ONE: Self;
}
}

That is the whole trait: identity constants and the accumulator type. f16 and bf16 accumulate in f32. i8 accumulates in i32. f32, f64, and the complex types accumulate in themselves. No arithmetic lives on Scalar itself. All real math happens vectorized in SimdOps, or in per-family scalar epilogues, so adding an element type never drags in a scalar arithmetic surface. The refinement traits Float, NarrowFloat, and ComplexFloat layer on the few extra capabilities specific paths need. What scalar.rs deliberately does not know: that SIMD exists. It has no idea its constants will end up broadcast into vector registers.

gemmkit/src/simd.rs forms the load-bearing wall. The simd/ backends join it: avx512.rs, fma.rs, neon.rs, scalar.rs, wasm.rs, and the complex glue in complex.rs. 3 traits split the job. Simd is a zero-sized ISA token. Examples include Avx512F, Fma, Neon, ScalarTok, Simd128, and the dot-capable Avx512Vnni and Avx512Bf16. Its sole method is vectorize, the #[target_feature] trampoline. It puts runtime-selected intrinsics into feature-enabled codegen. SimdOps<T> is the thick per-element-type vocabulary: register type, LANES, load/store/broadcast/mul/add/fma/reduce, and the overridable accumulate_tile schedule. KernelSimd<L, R, A, O> on top of them is the widen/narrow seam that makes mixed precision work without a driver branch. What this module deliberately does not know: anything above it. Its module doc states it depends only on scalar and core, so it could be split into its own crate unchanged. SimdOps has no idea what a micropanel, a cache, or a GEMM is.

L1: the mechanical copy, pack.rs

gemmkit/src/pack.rs holds the 2 shared packing primitives that turn a strided A or B region into contiguous, microkernel-sized panels. These are the copies the kernel families of L4 delegate their pack hooks to. The complex family’s plane-splitting pack is the one exception, and it lives with its family.

The 2 primitives are pack_panels and pack_kgroup_panels. pack_panels is the micropanel-major copy. LHS panels are mr rows tall, and RHS panels are nr columns wide. It is the same routine for both, with the “leading” and “depth” strides swapped. Tails are zero-filled, and a cache-blocked transpose handles strided sources. pack_kgroup_panels is the k-group-interleaved variant the dot-product families use.

What pack.rs deliberately does not know: where its output goes. The same routine fills a transient per-call scratch region, a shared parallel pack buffer, and a caller-held PackedRhs that lives for the whole process. pack.rs never sees a Workspace, a worker, or a lifetime, only dst, src, and strides. That indifference is what makes the prepacked path byte-identical to the per-call path. Depending only on scalar, it names no family and no cache, which is why it can sit this low.

L2: work distribution, parallel.rs

gemmkit/src/parallel.rs owns 3 things.

First, the Parallelism enum: Serial, or Rayon(n) with Rayon(0) meaning auto.

Second, workload-aware worker-count resolution. A serial gate applies below a total-work threshold. An explicit count is honored, but capped. An auto count scales with the total work m*n*k, rather than jumping straight to all cores. A separate bandwidth rule covers the memory-bound matrix-vector shapes.

Third, the demand-driven machinery. JobCursor is a lock-free atomic cursor workers pull contiguous chunks from. The job_grain and packed_block_grain knobs size those chunks. The for_each_worker fork-join is the barrier the higher layers use. parallel.rs also provides Ptr, the Send + Sync pointer shim that lets raw pointers cross into rayon closures.

Depending only on tuning, parallel.rs sits below everything that calls it. The same worker vocabulary serves kernel, driver, special, and dispatch alike. What it deliberately does not know: what a job is. JobCursor hands out index ranges over an abstract count. Nothing in this file mentions tiles, matrices, or families. That is why the same cursor later schedules driver tiles, B-pack panels, A-pack row blocks, and gemv row panels, without distinguishing between them. More in Parallel Execution.

L3: the machine model, cache.rs

gemmkit/src/cache.rs and its backends (cache/cpuid.rs, cache/sysfs.rs, cache/sysctl.rs) answer 2 questions. What does the cache hierarchy look like? What blocking follows from it? Detection is a best-effort fallback chain that cannot fail. It tries CPUID on x86, then Linux sysfs, then macOS sysctl, then a static default calibrated on a Zen5 part. #[cfg] only ever picks the sniffing method, never the values, and the result is memoized once in Machine.

blocking() then computes (MC, KC, NC) analytically from the BLIS model. KC is sized so the A and B micropanels coexist in L1. MC is sized so the A macro-panel fits L2. NC is sized so the B macro-panel fits L3. The key types are Level (with its carefully documented shared_by contention field), CacheTopology, and the blocking result. What this layer deliberately does not know: the thread count. blocking() has no worker parameter, and that omission is load-bearing. Thread-count-independent blocking is the mechanism behind the reproducibility contract described in Design Goals. Full detail in Blocking and the Cache Model.

L4: the operation-family seam, kernel.rs

gemmkit/src/kernel.rs and kernel/ (float.rs, mixed.rs, int.rs, complex.rs, epilogue.rs) define KernelFamily. This is the bundle of everything that distinguishes one kind of GEMM from another. A family bundles the Lhs, Rhs, Acc, and Out types. It bundles the pack layout too: pack_lhs and pack_rhs, which delegate to the L1 primitives. It bundles the microkernel, microkernel_epi. It also bundles constants like OUT_IS_ACC and DEPTH_MULTIPLE, which tell the driver how to block for the family.

FloatGemm<T> is the baseline. MixedGemm, IntGemm/IntGemmVnni, and ComplexGemm are siblings that reuse the driver unchanged. This layer also owns the Epilogue trait, with its zero-cost Identity. It owns the AlphaStatus and BetaStatus enums too. The driver precomputes both, so the microkernel never compares floats. What a family deliberately does not know: its own tile size. MR_REG and NR are const generics on the microkernel method, chosen per (type, ISA) at the dispatch site three layers up. The family compiles for any geometry, and a new tile is a new instantiation, never a new type.

L5: the engine, driver.rs

gemmkit/src/driver.rs is the one blocked loop nest that serves every family. It has the BLIS-order jc -> pc -> flat job list structure. It makes the adaptive packing decisions. It can pack B per depth slice, or read B in place. It can pack A per worker, pack A through a shared pre-pass, or not pack A at all. It also handles prepacked-RHS consumption, through the pack_rhs_full layout the prepack API reuses. That reuse is why prepacked and plain GEMM produce identical panel bytes.

Its public faces are run, run_epilogue, run_packed_rhs, and run_packed_rhs_epilogue, all funneling into the private run_inner. What it deliberately does not know: any concrete element type or ISA. The whole file is generic over Fam: KernelFamily and a KernelSimd token. It never names f32, never names AVX-512, and never branches on element type. That is the open/closed property. Adding a family or an ISA leaves this file untouched. gemmkit/tests/open_closed.rs proves this. It drives the driver with a second, trivial family the crate does not ship.

L6: the reroutes, special.rs

gemmkit/src/special.rs and special/ (gemv.rs, small_k.rs, small_mn.rs, batched.rs) hold the paths for shapes the register-tiling driver fits poorly. These are matrix-times-vector, low-depth GEMM, small-m,n long-k inner products, and the batched orchestration layer. All sit behind the same public entries and are covered in Special Paths. What a special path deliberately does not know: why it was chosen. The gates, gemv_threshold, small_k_threshold, and small_mn_dim, live in the dispatch layer above and the tuning module beside it. small_k::run cannot even tell whether it is serving gemm, gemm_fused, or gemm_map, because the epilogue arrives as an opaque generic parameter. The batched path is the one place a layer reaches back up. batched.rs forwards each element through dispatch::execute (L7), so it inherits the whole ladder above. That is the single annotated exception discussed at the end of this page.

L7: runtime ISA selection, dispatch.rs

gemmkit/src/dispatch.rs and dispatch/ (isa.rs, float.rs, mixed.rs, int.rs, complex.rs) turn “which kernel should this machine run” into a one-time decision. Each element type has one OnceLock<Dispatched<T>> slot. Feature detection runs once. The winning monomorphized entry points (plain, prepacked, fused) get cached, along with the tile geometry. Every later call is a plain indirect call through a typed function pointer, with no transmute and no AtomicPtr<()>.

This layer also owns the Task<T> problem descriptor, the degenerate-case handling in execute, the orientation normalization orient_transpose, and the special-path gates. It also owns the GEMMKIT_REQUIRE_ISA pin, which forces, or fails loudly on, a specific kernel. What it deliberately does not know: where the pointers in a Task came from. Checked slice views and raw unchecked pointers arrive identical. Validation happened above or not at all, and dispatch neither knows nor cares. See SIMD Tokens and ISA Dispatch and the user-facing Runtime ISA Dispatch.

L8a: the public boundary, api.rs

gemmkit/src/api.rs and api/ (batched.rs, cplx.rs, fused.rs, int8.rs, map.rs, packed.rs) define several things. These are the MatRef/MatMut strided views, the per-family safe entries, and the validate_gemm_views panic catalog. The safe entries come in *_with (caller-owned workspace) and *_unchecked (raw engine) variants. This layer also lowers views into Tasks.

What it deliberately does not know: everything below dispatch. The API layer cannot see which ISA will run, what blocking will be chosen, or whether packing will happen. After validation it hands a Task to dispatch::execute and its job is done. Symmetrically, MatRef never appears below this layer. The rest of the crate speaks only pointers and strides.

Why the arrows only point down

The dependency direction is the architecture’s one hard rule. Each layer is driven by the layers above it, and knows nothing about them. simd depends only on scalar and core. The driver never names an element type or ISA. Nothing below L7 knows dispatch exists. Nothing below L8a has ever heard of a slice.

The rule has exactly one deliberate, annotated exception. special/batched.rs (L6) forwards each batch element back through dispatch::execute (L7). This re-entry lets every element inherit the same driver, small-k, small-mn, and gemv routing a standalone gemm call would take. It avoids a second dispatch ladder maintained by hand. It is an upward arrow by design, and the only one in the crate.

3 payoffs justify the discipline. First, extension cost. Because knowledge only flows downward, a new ISA, element type, or family plugs in at its own layer. Everything below it is provably untouched. The seams described in Extension Points work only because no lower layer could have special-cased what sits above it. Second, review locality: to audit the microkernel, read kernel/float.rs and the SimdOps contract, nothing else. To audit the scheduling, read driver.rs and parallel.rs. Third, testability: lower layers are exercised in isolation. SIMD conformance tests check every token against scalar models, and the open/closed test drives the driver with a foreign family. This is what makes the correctness story in Testing and Verification tractable.

The 2 cross-cutting modules

2 modules sit beside the stack rather than in it, because every layer needs them and neither depends on anything above core/alloc.

gemmkit/src/tuning.rs is the unified knob surface. Every heuristic threshold in the engine lives here: the serial/parallel gate, pack gates and strides, the special-path thresholds, scheduler grains, and blocking caps. Each one resolves in this order: per-call argument, then programmatic setter (tuning::set_*), then environment variable (GEMMKIT_*), then compiled default. Env vars are read once and cached. A malformed value warns on stderr and falls back rather than panics, because a perf-knob typo must not crash the process. The full set of GEMMKIT_* names is enumerated in the tuning::knob_env_names registry. The out-of-crate consumers, the gemmkit-tune sweep table, the knob property tests, and the fuzz setters, assert their lists against it. So a new knob cannot silently escape coverage. The user-facing tour is Tuning Knobs.

gemmkit/src/workspace.rs is the scratch-memory story. Workspace is a growable 64-byte-aligned buffer. Workspace::regions carves it into per-worker (or per-row-block) LHS regions plus one shared RHS region, with fail-closed overflow checks at the element-to-byte chokepoint. Under std a re-entrancy-safe thread-local pool supplies the default, so plain gemm allocates at most once per thread. The *_with entries thread a caller-owned workspace through instead, giving zero heap allocation after the first sufficiently large call. Without std, each call uses a fresh workspace. Details in Packing and Workspaces.

Life of a GEMM Call

The previous page described the stack at rest. This page follows a single call through it. The specimen is the quick-start example from the crate docs:

#![allow(unused)]
fn main() {
use gemmkit::{gemm, MatRef, MatMut, Parallelism};

// 2x3 * 3x2 = 2x2, all row-major
let a = [1.0_f32, 2.0, 3.0, 4.0, 5.0, 6.0];
let b = [7.0_f32, 8.0, 9.0, 10.0, 11.0, 12.0];
let mut c = [0.0_f32; 4];
gemm(
    1.0,
    MatRef::from_row_major(&a, 2, 3),
    MatRef::from_row_major(&b, 3, 2),
    0.0,
    MatMut::from_row_major(&mut c, 2, 2),
    Parallelism::Serial,
);
assert_eq!(c, [58.0, 64.0, 139.0, 154.0]);
}

This toy shape takes one of the early exits below. The walk keeps 2 problems in mind: the 2x2x3 example above, and a 2048x2048x2048 f32 product on an AVX-512 machine. That larger product goes all the way down, through every layer. Here is the route, compressed:

gemm(alpha, A, B, beta, C, par)
  |  validate_gemm_views: shapes, bounds, aliasing     [api.rs]
  v
Task<T>: raw pointers + isize strides
  |  m == 0 || n == 0        -> return                 [dispatch.rs]
  |  k == 0 || alpha == 0    -> C <- beta*C, done
  v
memoized per-type kernel (OnceLock fn pointer)
  |  gemv shape (m==1||n==1) -> special/gemv.rs        [dispatch/float.rs]
  |  orient: row-major-ish C -> compute C^T = B^T*A^T
  |  small m,n + long k      -> special/small_mn.rs
  |  k <= small_k_threshold  -> special/small_k.rs
  v
driver::run                                            [driver.rs]
  jc over NC -> pc over KC (never parallel)
    -> flat job list (ic row-block x jt column-tile),
       workers drain a shared JobCursor, pack A/B adaptively
  v
Fam::microkernel_epi: MR x NR tile in registers        [kernel/float.rs]
  alpha/beta epilogue store (vector fast path | scratch drain)

Stage 1: validation and lowering

gemm itself is one line. It borrows the thread-local workspace and forwards to gemm_with (gemmkit/src/api.rs). gemm_with runs validate_gemm_views, the full panic catalog from Design Goals. Shapes must agree. Every view must stay inside its slice. C must address each (i, j) uniquely. C must overlap neither input.

Then the views dissolve. Everything below this point speaks Task<T>: a Copy struct of m, k, n, alpha/beta, and 3 raw pointers with isize row/column strides. Transposition never exists as a flag. A transposed view is just swapped strides. When beta == 0, the contract says C is never read, so it may arrive uninitialized.

The unsafe boundary is crossed exactly here, justified by the validation that just ran. gemm_unchecked enters one step later, with the caller carrying that justification instead.

Stage 2: dispatch early exits

dispatch::execute (gemmkit/src/dispatch.rs) handles the degenerate algebra while the element type is still concrete:

#![allow(unused)]
fn main() {
if task.m == 0 || task.n == 0 {
    return;
}
// k == 0 or alpha == 0 => the A*B term vanishes: C <- beta*C only
if task.k == 0 || task.alpha == T::ZERO {
    T::scale_c(task.beta, task.c, task.m, task.n, task.rsc, task.csc);
    return;
}
T::dispatch(task, par, ws);
}

An empty output means nothing to do. A vanished A*B term (k == 0 or alpha == 0) degrades the call to a C <- beta*C scale. That scale never reads A or B. Within it, beta == 0 stores zeros without reading C, which keeps the uninitialized-C contract honest.

Only a real product reaches T::dispatch, which reads the per-type OnceLock slot. On first use, the selection ladder probes CPU features, honoring a GEMMKIT_REQUIRE_ISA pin, which panics rather than falls back. It then caches the winning monomorphized entry points plus the tile geometry. Every later call is a single indirect call. On the AVX-512 machine, f32 resolves to run_typed::<f32, Avx512F, 2, 12>, the 32x12 tile.

Stage 3: routing in run_typed

run_typed (gemmkit/src/dispatch/float.rs) is a short gauntlet of gates, each rerouting a shape the register-tiling driver would serve poorly.

First, gemv. If m == 1 || n == 1, and the path is not capped off through GEMMKIT_GEMV_THRESHOLD, the call goes straight to special/gemv.rs. This happens before orientation normalization, in the user’s original frame. gemv resolves its own orientation. It treats the m == 1 case as the transposed rows x k problem, and it partitions output rows itself.

Everything else is orientation-normalized by orient_transpose. If C is row-major-ish (|csc| < |rsc|), the dispatcher rewrites the problem as its transpose: C^T = B^T * A^T. This swaps m with n, the A/B pointers and strides, and rsc with csc.

The identity is free. No data moves, only the descriptor changes. It buys a strong invariant. After this point, the output’s row stride is the small one (rsc == 1 for a fully contiguous C). Each output column then occupies consecutive memory, and the kernel walks down contiguous columns.

The microkernel’s fast store path needs exactly that: rsc == 1, so it can use vector stores of LANES consecutive rows in a column. Every layer below optimizes for one orientation instead of 2. The all-row-major 2048-cube hits this swap. The engine actually computes C^T, and nobody below dispatch knows.

Then 2 more gates apply to the normalized task.

A small-m,n shape goes to special/small_mn.rs. Both dimensions must be at or below small_mn_dim, with a contraction longer than small_k_threshold. There, each output element is one horizontal SIMD dot. This is zero-copy when both operands stream unit-stride along k. When one operand is strided (k > small_mn_pack_min_k), it goes through a pack tier that copies just the offending operand.

A small-k shape, k <= small_k_threshold (16 on x86, 8 on aarch64 by default), goes to special/small_k.rs. This computes the whole product as one in-place depth panel over the microkernel, with no blocking or packing setup.

Whatever passes all the gates, and the 2048-cube does, enters driver::run. The driver states its preconditions: m, n, k > 0, alpha != 0, and orientation normalized.

Stage 4: the driver loop nest

driver::run forwards to run_inner (gemmkit/src/driver.rs) with the zero-cost Identity epilogue. The fused entries land in the same function, with a real epilogue instead. The driver is generic over the family and the ISA token. For this call, that is FloatGemm<f32> and Avx512F, with mr = MR_REG * LANES = 32 and nr = 12.

Blocking comes first. cache::topology().blocking(mr, nr, sizeof_lhs, m, n, k) yields (MC, KC, NC) from the BLIS cache model. These are sized in packed-input elements, sizeof(Lhs), not the accumulator, so narrow types get deeper blocks. The loop nest then runs in BLIS order:

• jc over NC: column blocks, sized so the packed B macro-panel stays L3-resident.
• pc over KC: depth slices. This loop is never parallel. All depth slices accumulate into the same C tiles. Parallelizing depth would mean a synchronized read-modify-write on C, or split reductions. Keeping depth serial lets every output element be reduced start to finish by one worker, which is half of the reproducibility contract. beta participates only on the first slice (pc == 0). Later slices run with an effective beta of one, and accumulate. For mixed-precision families (OUT_IS_ACC = false), there is exactly one slice, kc = k, so the running sum never rounds through the narrow output type.
• A flat 1-D job list: inside each depth slice, the remaining work is n_mc row blocks times n_nt column tiles. These flatten into n_jobs = n_mc * n_nt indices. Workers pull contiguous chunks from a shared, lock-free JobCursor on demand. There is no static partition, so a faster core absorbs proportionally more work. The chunk grain oversamples the worker count (job_grain). The packed-LHS path uses a row-block-aligned packed_block_grain instead, so chunks never straddle a pack boundary. The worker count itself comes from par.resolve(m*n*k, n_jobs). This is work-based: it scales with the total work m*n*k over a per-worker floor, rather than jumping straight to all cores. If that count would leave the job list shallower than a few chunks per worker, the driver first shrinks mc. This only cuts more, smaller row blocks, so it cannot move a result bit. Shrinking mc deepens the list before the cursor hands work out.

Packing is adaptive, and each side decides separately.

B is packed once per depth slice, when m clears rhs_pack_threshold. The packed panel is reused across all n_mc row blocks, so the copy only pays off when that reuse is high. Otherwise B is read in place, through its original strides. When packing does happen, it is itself parallel. Workers pull nr-wide column panels from a cursor. The for_each_worker join doubles as the write-before-read barrier. Packed B is the one buffer all compute workers share, so this barrier matters.

A has 3 modes. Each worker can pack the row block it is working on into its own private workspace region. This is forced when rsa != 1, or when the block is a partial mr multiple. It is chosen otherwise when per-worker column reuse, or a TLB-hostile column stride, makes it pay off. On large parallel problems, a shared pre-pass can instead pack each row block exactly once. It packs into a per-block region behind its own barrier (the shared_lhs_mnk gate). This eliminates redundant per-worker packing. When reuse is too low to pay for any copy, A is read in place.

Sizing for these regions happens up front, through Workspace::regions, with the fail-closed overflow checks noted earlier. On a no-pack route, the workspace is not even touched.

For each job, the worker resolves its A panel, packed or in-place. It locates the B panel: per-call packed, prepacked buffer, or in-place. It then calls the microkernel for every mr-row strip of the block. All of this runs inside simd.vectorize, so the entire strip executes in target-feature codegen.

Stage 5: the microkernel and its store

Fam::microkernel_epi (gemmkit/src/kernel/float.rs, microkernel_impl) computes one MR x NR tile. For this call, that is 32x12 f32 values, held in 24 ZMM accumulator registers as a [[Reg; MR_REG]; NR] array.

A full-width tile runs SimdOps::accumulate_tile, the ascending-k fused-multiply-add schedule. This is the seam a load-bound ISA like NEON overrides, with a software-pipelined variant that reorders loads, never arithmetic. An edge column tile instead takes a runtime-bounded loop. That loop reads exactly nr_eff columns, so an unpacked B is never read past its last real column.

Then alpha folds into the accumulators. This step is skipped entirely when alpha == 1, thanks to the AlphaStatus the driver precomputed.

The store is where beta and the epilogue live. It has 2 routes.

The fast path fires for a full tile with unit output row stride: mr_eff == mr && nr_eff == NR && rsc == 1. The orientation normalization from stage 3 is what makes this common. Each accumulator register combines with C directly. It is stored as-is for beta == 0 (C unread), added for beta == 1, or fused-multiply-added for a general beta. The result is written back with vector stores.

Edge tiles and strided outputs take the general path instead. All accumulators drain into a stack scratch tile, a SCRATCH_LEN array in the worker’s frame, with no allocation. A scalar loop then applies the same beta arithmetic element-wise, through whatever strides C has.

Plain gemm threads the Identity epilogue through all of this. Every epilogue hook is gated on !E::IS_IDENTITY, an associated const. The guards fold away at monomorphization, so the emitted kernel is byte-for-byte the pre-epilogue code.

A fused call, such as gemm_fused, gemm_map, or requantize, runs the same engine with a real epilogue. That epilogue fires only when last_k is true, on the final depth slice, once per output element. That story continues in Epilogue Fusion.

The short way home

The 2x2x3 example never saw most of this. It entered execute with m, n, k all positive and alpha == 1. It reached run_typed, and failed the gemv gate (n != 1, m != 1). It was orientation-swapped, then failed the small-m,n gate, since k = 3 is not a long contraction. With k = 3 <= 16, it took special/small_k.rs. That route is one in-place depth panel over the same microkernel, with no blocking, no packing, and no workspace traffic.

The 2048-cube took the full driver, with parallel B-packing. At Parallelism::Rayon(0), its worker count scaled to its total work.

Same entry point, same result contract, 2 very different journeys. The layers below decide, and the caller never has to. For the deeper mechanics of each stage, see Blocking and the Cache Model, Packing and Workspaces, Parallel Execution, and Special Paths.

SIMD Tokens and ISA Dispatch

gemmkit picks its instruction set at runtime. That decision collides with how Rust compiles SIMD intrinsics. AVX and AVX-512 intrinsics compile correctly only inside a context where the target feature is enabled. Normally that context comes from a #[target_feature(enable = "...")] attribute on the enclosing function. The program only knows which features are safe to enable once it runs on a concrete CPU. The microkernel is one generic function shared by every instruction set, so no single attribute can sit on it.

This page explains 2 things. First, how the L0 SIMD layer (gemmkit/src/simd.rs and gemmkit/src/simd/) resolves that tension with zero-sized ISA tokens and a trampoline function. Second, how the L7 dispatch layer (gemmkit/src/dispatch.rs and gemmkit/src/dispatch/) selects and caches the winning kernel.

ISA tokens and the vectorize trampoline

An ISA token is a zero-sized type that stands for one instruction-set choice. On x86 the tokens are Fma (AVX2 + FMA) and Avx512F, plus the dot-kernel variants Avx512Vnni and Avx512Bf16. On aarch64 the token is Neon. On wasm32 the token is Simd128. ScalarTok exists on every platform as the portable floor.

Each token implements the Simd trait. Its only method is vectorize, which runs a closure with this token’s target features enabled. The code below shows the entire mechanism, from gemmkit/src/simd/fma.rs.

#![allow(unused)]
fn main() {
/// AVX2 + FMA ISA token
#[derive(Copy, Clone, Default)]
pub struct Fma;

impl Simd for Fma {
    #[inline(always)]
    unsafe fn vectorize<R>(self, f: impl FnOnce() -> R) -> R {
        #[target_feature(enable = "avx2,fma,f16c")]
        unsafe fn inner<R>(f: impl FnOnce() -> R) -> R {
            f()
        }
        // SAFETY: the caller of `vectorize` (the runtime dispatcher) guarantees
        // the CPU supports avx2+fma(+f16c); `inner` then establishes the codegen
        // context, and `f` inlines into it
        unsafe { inner(f) }
    }
}
}

The trick is the direction of inlining. inner is a tiny function with a #[target_feature] attribute. The closure f inlines into it. f holds the packing loops and the microkernel calls, and every one of those is built from #[inline(always)] primitives. Every intrinsic therefore lands in a codegen context where the feature is enabled. No attribute ever touches the generic kernel itself.

The unsafe contract has exactly one obligation. The caller must guarantee the CPU really supports the token’s features. The runtime dispatcher establishes this once per process.

This is the same pattern pulp and faer use. It works the same way for the serial path and for rayon worker closures. The driver wraps each column strip of microkernel calls in simd.vectorize(|| ...). This amortizes the trampoline overhead over many tiles.

ScalarTok’s vectorize is just f(). It has nothing to enable. This is what makes the scalar path run everywhere, including under Miri.

SimdOps: the per-type vocabulary

L0 builds on 3 traits, not 2. Simd is the ISA token trait above. SimdOps<T> is the per-element-type vocabulary this section covers. KernelSimd<L, R, A, O> is a third trait. It widens loads and narrows stores when a family’s input type, accumulator type, and output type are not all the same. This section covers it near the end.

The token itself knows nothing about element types. All the raw operations live on SimdOps<T>, implemented once per (ISA, T) pair. It names the register type Reg, the lane count LANES, and every primitive the microkernel needs. The token and the element type are decoupled, so LANES varies with the pair. f32 gets 8 lanes under Fma and 16 lanes under Avx512F. f64 gets half as many lanes as f32 on the same token.

The vocabulary is deliberately thick. The basic operations are zero, splat, loadu, storeu, mul, add, and the fused mul_add. Its subtractive partner is fnma, which computes c - a*b. The complex kernel needs fnma for one of its accumulation terms. The vocabulary also has the horizontal reduce_sum, used by the gemv and dot-product epilogues.

On top of those sit a few more primitives. max and min exist only on the real-float tokens, for the fused ReLU and clip epilogues. The LANE_FMA flag and its fma_bvec method give NEON a lane-indexed FMA path. It loads a block of RHS columns as one vector, instead of issuing a separate splat per column. accumulate_tile is the GEMM inner loop itself. Its portable default schedule already compiles down to the canonical register-blocked kernel on any out-of-order core.

The complex split kernel has its own seam here too, cplx_microkernel. The dot kernels have theirs as well, dot_accumulate, which lives on KernelSimd rather than on SimdOps. KernelSimd<L, R, A, O> is the seam a family drives on when its input types, accumulator type, and output type are not all equal. One example is an f16 input with an f32 accumulator. It widens a narrow input load into the accumulator type and narrows an accumulator value back down on store. A homogeneous family has all 4 types equal. It gets a KernelSimd implementation for free through a blanket impl. It never needs any per-ISA code of its own. See Scalars and Kernel Families and Dot Kernels and the Deep-K Twin for more on both seams.

The thickness is the point. matrixmultiply’s per-ISA trait is thin, so each instruction set has to reimplement the kernel from scratch. Here every primitive the kernel needs sits behind SimdOps, so the microkernel is one generic function over every ISA. Adding an instruction set costs a new token, its SimdOps impls, and one line in each dispatch ladder. The simd module depends only on crate::scalar and core. It has no reverse dependency on the kernel, driver, or cache layers, so the whole abstraction could move into its own crate unchanged.

The dispatch layer

Dispatch turns the question “which token applies” into a one-time decision. Each dispatched element type owns one OnceLock slot holding a Dispatched<T> descriptor. f32 and f64 use gemmkit/src/dispatch/float.rs. f16 and bf16 use dispatch/mixed.rs. i8 uses dispatch/int.rs, with its own IntDispatched and IntRequantDispatched shapes, because those types are heterogeneous. c32 and c64 use dispatch/complex.rs. The code below, lightly trimmed, is from dispatch/float.rs.

#![allow(unused)]
fn main() {
#[derive(Copy, Clone)]
pub(super) struct Dispatched<T> {
    pub(super) run: GemmFn<T>,
    pub(super) run_packed: PackedFn<T>,
    #[cfg(feature = "epilogue")]
    pub(super) run_fused: FusedFn<T>,
    #[cfg(feature = "epilogue")]
    pub(super) run_packed_fused: PackedFusedFn<T>,
    pub(super) mr: usize,
    pub(super) nr: usize,
    pub(super) depth_multiple: usize,
}
}

The slot caches the winning monomorphized entry points: the plain kernel, the prepacked-RHS kernel, and, under the epilogue feature, their fused twins. It also caches the microtile geometry (mr, nr) and the family’s depth_multiple. The geometry is cached so prepack_rhs can size a buffer through the same ISA choice the consuming call will make. depth_multiple lets the bf16 prepack path round its packed depth to match the dot kernel’s layout. Everything here is a typed function pointer. There is no transmute and no AtomicPtr<()>.

A call flows through a fixed chain. gemm calls dispatch::execute, which handles the degenerate cases. dispatch::execute calls T::dispatch, which reads the memoized slot. The slot resolves to one indirect call into a wrapper, such as gemm_f32_avx512f. That wrapper instantiates the shared generic entry as run_typed::<f32, Avx512F, 2, 12>.

Selection runs once, inside the OnceLock initializer. It first honors any GEMMKIT_REQUIRE_ISA pin, covered below. After that, the auto ladder on x86 probes avx512f first, then avx2 plus fma, and falls back to scalar last. On aarch64, NEON is the baseline. The architecture makes NEON mandatory, so no probe is needed there.

On wasm32 there is no runtime feature detection at all. simd128 is chosen at compile time, through cfg(target_feature = "simd128"). The build must pass -C target-feature=+simd128, or it gets the scalar kernel instead. Scalar is the floor on every architecture.

Each per-type ladder adds its own gate on the same skeleton. The f16 FMA arm also requires f16c, for the vcvtph2ps and vcvtps2ph conversions. The bf16 ladder tries the avx512bf16 dot kernel before plain AVX-512F. The i8 ladder tries avx512vnni (with avx512bw) before the widen kernel.

2 build-mode details are worth knowing. Under std, feature detection uses is_x86_feature_detected!, and the result is memoized in the OnceLock. Without std, there is no runtime CPU detection, because raw-cpuid is gated on std. The probe macro degrades to cfg!(target_feature = ...), GEMMKIT_REQUIRE_ISA parsing degrades to Auto, and the select function runs on every call. Every branch inside it is now a compile-time constant, so it folds down to a direct choice. A no_std build simply runs whatever its compile-time target features guarantee. See no_std and WebAssembly.

Tile geometry as const generics

Besides the instruction encoding, the one thing that genuinely varies per (type, ISA) is the microtile shape. It is expressed as a pair of const generics, (MR_REG, NR), chosen at the dispatch site. It is never a new type, trait, or macro. MR_REG is how many registers tall the tile is, so the row count is MR = MR_REG * LANES. The table below covers f32.

f64 halves the lane count. The same (MR_REG, NR) pairs then yield 16x12 on AVX-512F, 8x6 on FMA, 8x4 on NEON, and 4x4 on simd128. The budgets are not accidents. NEON deliberately leaves about 11 registers free, as rename headroom for a wide out-of-order core to overlap loads with FMAs. simd128 stays at 11 live vectors, because LLVM’s wasm backend starts spilling past roughly 16. These comments live next to the wrappers in dispatch/float.rs, so the table above is the code, not an aspiration.

Pinning with GEMMKIT_REQUIRE_ISA

By default the best available ISA wins. Setting the environment variable GEMMKIT_REQUIRE_ISA forces exactly one kernel instead of the automatic choice. It accepts these values, case-insensitive:

• scalar
• fma (alias avx2)
• avx512f
• avx512vnni (alias vnni)
• avx512bf16 (alias bf16)
• neon
• simd128 (alias wasm)
• auto

Unset or empty means auto. An unrecognized value causes a hard panic, so a typo in a CI configuration cannot silently select the wrong thing. avx512vnni pins the i8 vpdpbusd dot kernel. avx512bf16 pins the bf16 vdpbf16ps dot kernel. For every other element type, both resolve to the plain AVX-512F path.

The defining behavior is that a pin never falls back. Selection panics if the CPU, or an emulator such as Intel SDE, does not report the required feature. It also panics if the requested ISA does not exist on the target architecture at all, such as neon off aarch64 or fma/avx512* off x86. The panic message names the missing feature. The rationale is CI honesty. A job that means to exercise a given kernel must fail loudly, rather than silently test a different one.

The simd128 pin earns its keep the same way on wasm. The target feature is an easily forgotten compile-time flag. Pinning turns a dropped flag from a silent scalar fallback into a build that refuses to run.

The value is read once. Selection is memoized in the per-type OnceLock, so the variable must be set in the process environment before the first GEMM call. Changing it afterward has no effect for the life of the process. See Runtime ISA Dispatch for the user-facing side of pinning, including CI recipes and how it interacts with the tuning knobs.

Scalars and Kernel Families

gemmkit multiplies several element types: f32, f64, f16, bf16, i8, Complex<f32>, and Complex<f64>. u8 also appears, but only as a requantize output. Every one of these types flows through the same driver, the same packing framework, the same cache model, and the same parallel scheduler.

2 traits carry all the variation. Scalar, at L0 (gemmkit/src/scalar.rs), answers what a type is and what it accumulates in. KernelFamily, at L4 (gemmkit/src/kernel.rs), answers what makes this kind of GEMM different from the others. The driver is generic over the family. It never branches on element type. This page walks through that split, and why it falls where it does.

Scalar: constants and an accumulator, nothing else

Scalar is deliberately tiny. The whole trait, from gemmkit/src/scalar.rs:

#![allow(unused)]
fn main() {
pub trait Scalar: Copy + Send + Sync + PartialEq + 'static {
    /// The type in which products are accumulated. `Self` for `f32`/`f64`
    type Acc: Scalar<Acc = Self::Acc>;
    /// The additive identity
    const ZERO: Self;
    /// The multiplicative identity
    const ONE: Self;
}
}

Scalar has no Add, no Mul, and no conversions. It has only the identity constants and the associated accumulator type. This omission is deliberate. All vectorized arithmetic lives in SimdOps (see SIMD Tokens and ISA Dispatch). The scalar arithmetic that an epilogue needs lives instead on narrow side traits, and only the types that need them implement those traits. Float covers f32 and f64, and, through num-complex’s operators, the complex types too. NarrowFloat covers the f16/bf16 widen and narrow conversions. ComplexFloat covers the real and imaginary accessors of the split complex kernel.

If Scalar itself carried arithmetic, every new element type would owe a full set of operations it may not actually have. i8 is the clearest case. It needs no arithmetic trait at all, because its kernel does everything through the SIMD seam and exact i32 integer operations.

Acc is the mixed-precision seam. The table is short.

The recursive bound Acc: Scalar<Acc = Self::Acc> pins the chain after one step: f16 -> f32 -> f32 -> .... Generic code can then name the accumulator’s accumulator without caring how narrow the input was. For the homogeneous types, the Acc = Self branch collapses at compile time and costs nothing.

KernelFamily: everything that distinguishes one GEMM from another

KernelFamily bundles the rest. It carries the 4 element types (Lhs, Rhs, Acc, Out), the pack layout (pack_lhs/pack_rhs, which write micropanel-major panels), and the microkernel.

3 associated constants shape how the driver treats a family. OUT_IS_ACC says whether a running partial sum may round-trip through C between depth panels. This is the pivotal constant, covered below. FORCE_PACK_LHS and FORCE_PACK_RHS are set when packing performs a transform the kernel depends on, such as complex conjugation or dot-kernel interleaving. In that case, the driver must never read the operand in place. DEPTH_MULTIPLE is the instruction-group depth of a dot kernel. It is 1 for every other family.

A family overrides exactly one of 2 microkernel methods. A non-fusing family overrides the plain microkernel and inherits the default microkernel_epi. That default forwards to microkernel, after it asserts E::IS_IDENTITY, a fail-closed guard. A real epilogue reaching a family that cannot fuse panics, instead of being silently dropped. A fusing family, such as the float, mixed, and requantizing families, overrides microkernel_epi instead. It threads the epilogue through its own store, and its plain microkernel method is then dead code, keeping the default unreachable! body.

Tile geometry is deliberately not on the trait. (MR_REG, NR) is a pair of const generics, chosen per (family, ISA) at the dispatch site. A new tile is therefore a new instantiation of that pair, not a new type.

The payoff shows up in the driver’s signature. driver::run::<Fam, S, MR_REG, NR> is generic over the family and the ISA token. It calls Fam::pack_lhs, Fam::pack_rhs, and Fam::microkernel_epi, and it contains not one if on an element type. Adding a kind of GEMM means writing a new family. It never means touching the driver.

The family roster

10 family types ship today. They fall into generations: homogeneous, widen, dot, requantize, and complex. Reading them in that order makes the seams visible.

FloatGemm is the reference point. It is homogeneous, with one generic microkernel_impl shared by every ISA and every tile.

The mixed and integer families introduce Acc != Lhs. They lean entirely on the widen/narrow seam covered below. The dot families, Bf16DotGemm and IntGemmVnni, go further. Each swaps in an interleaved pack layout and a hardware dot instruction. The f32-output twins exist so a deep contraction can re-block. Dot Kernels and the Deep-K Twin covers all of this.

The requantizing variants bolt an exact i32 -> i8/u8 requantize onto the integer accumulation. They are part of the fusion story in Epilogue Fusion.

ComplexGemm keeps Acc = T, so complex alpha/beta thread through the driver unchanged. Its hot loop instead runs on the real component, through a dedicated seam. The Complex Split Kernel covers that seam.

This page stays at the roster level. The deep dives live in those other pages.

KernelSimd: the widen/narrow seam

The driver’s bound on the ISA token is S: KernelSimd<Fam::Lhs, Fam::Rhs, Fam::Acc, Fam::Out> (gemmkit/src/simd.rs). KernelSimd<L, R, A, O> extends SimdOps<A>, so it accumulates in A. It adds 4 moves a family needs at the type boundary.

load_lhs loads LANES LHS values and widens them to an A register. splat_rhs widens one RHS scalar and broadcasts it. load_out widens output values for the beta != 0 read of C. store_out narrows an A register to LANES output values, and rounds to nearest-even when it actually narrows.

The homogeneous case costs nothing. A blanket impl, KernelSimd<A, A, A, A> for S: SimdOps<A>, forwards all 4 methods to plain loadu, splat, and storeu. So FloatGemm<f32> and its relatives need zero per-ISA code.

A mixed family instead adds per-ISA impls. Its loads genuinely widen, such as f16 -> f32 through vcvtph2ps, or i8 -> i32 through sign extension. Its store_out genuinely narrows. Coherence comes free here. The all-equal blanket and a mixed impl with L != A can never describe the same types.

2 further impl groups are derived, rather than hand-written per ISA. The requant blankets cover Out = i8 or u8. They forward the accumulate side to the <i8, i8, i32, i32> impl. The f32-output twins cover <N, N, f32, f32> for N = f16 or bf16. Those are written as 2 concrete heads, rather than one generic blanket over N. A generic blanket could not rule out colliding with the homogeneous blanket at N = f32.

KernelSimd also hosts 2 more seams. dot_accumulate is the dot seam. Only dot-capable tokens override it, and its default is unreachable!. requant_store is the vectorized requantize store, following the same pattern.

The constant that ties this seam to the driver’s blocking is OUT_IS_ACC. The driver normally accumulates across k by splitting it into kc panels. The partial sum round-trips through C, with beta = 1 after the first panel. That round-trip is exact only when Out == Acc.

When the output is narrower than the accumulator, that round-trip would round to 16 bits at every panel boundary. So a narrow family declares OUT_IS_ACC = false. The driver then responds with kc = k, one depth panel that spans the entire contraction. The whole contraction accumulates in f32 registers, and the result rounds to the narrow output exactly once, at the end.

That single-rounding guarantee is what makes the mixed-precision results defensible. It has one cost. A single panel means its RHS micropanel can outgrow the L2 cache when k is large. The f32-output twins exist to pay that cost down. Dot Kernels and the Deep-K Twin covers how they do it.

The open/closed proof

The claim that the family seam is open for extension is not just prose. gemmkit/tests/open_closed.rs enforces it. This is an integration test that lives outside the crate, so it sees only the public API.

The test declares NaiveFloat, a deliberately naive second float family that shares nothing with FloatGemm. It reimplements micropanel packing from scratch, because the crate’s internal pack helper is not visible to it. This is exactly the situation a third party would be in. NaiveFloat also supplies a plain scalar triple-loop microkernel.

The test then drives the unchanged generic driver, driver::run::<NaiveFloat, ScalarTok, 4, 4>, on a 40x33x28 problem. It checks the result against an f64 reference.

The test’s main value is that it compiles at all. A second family drove the driver with no edit to driver.rs or pack.rs. It used nothing but public items: gemmkit::kernel::KernelFamily, gemmkit::simd::ScalarTok, gemmkit::driver::run, Workspace, and Parallelism.

Any refactor that closes the seam breaks this file first. Examples include a driver branch on a concrete family, a newly required private helper, or a leaked internal type in the trait’s signature. Such a refactor breaks this test before it breaks a downstream user.

Testing and Verification covers the wider testing story, including how the real families are cross-checked against oracles. Extension Points covers what third parties can build on this seam.

Dot Kernels and the Deep-K Twin

Most kernel families consume the contraction one depth step at a time. Each step loads a column of the packed A panel, and broadcasts one element of packed B. It then issues one FMA, or widening multiply-add, per accumulator register.

2 AVX-512 extensions break that rhythm. Each folds several depth steps into a single instruction. VNNI’s vpdpbusd multiplies 4 consecutive i8 depth steps into each of 16 i32 lanes. AVX-512 BF16’s vdpbf16ps folds 2 consecutive bf16 depth steps into each f32 lane.

An instruction that consumes several depth steps at once wants those steps adjacent in memory. For floats, it also changes how the accumulation rounds. So gemmkit gives dot kernels their own kernel families and ISA tokens, instead of hiding them behind a branch in the shared microkernel.

This page walks through 3 things. First, the 2 seams that carry dot kernels. Second, the 2 concrete dot kernels themselves. Third, the deep-contraction route that both narrow families share.

Why dot instructions get their own families

A dot kernel differs from its widen sibling in exactly 2 places. Each difference lands on a different extension axis of the engine. See Scalars and Kernel Families for the family and token split.

The pack layout is a family concern. KernelFamily::pack_lhs and pack_rhs take no ISA parameter, so a different interleave must key off the family instead. This is why Bf16DotGemm is a sibling of MixedGemm<bf16>, rather than a branch inside it.

The inner loop is a token concern. Only a CPU that actually has vpdpbusd or vdpbf16ps can run it. So the instruction lives behind a KernelSimd method, and only dot-capable tokens override that method.

The concrete dot families are IntGemmVnni, and its requantizing variant IntGemmVnniQ, with 4 depth steps per instruction. Bf16DotGemm has 2 depth steps per instruction. 2 f32-output twins, covered at the end of this page, round out the set.

DEPTH_MULTIPLE and the k-group pack

A family that folds Q depth steps per instruction declares const DEPTH_MULTIPLE: usize = Q. The default is 1. The contract, spelled out in gemmkit/src/kernel.rs, works as follows. The family’s packers write panels of width * kc.next_multiple_of(Q) elements, and pad the depth tail. The driver strides packed panels by that same padded depth, so both sides stay in lockstep. For every ordinary family, DEPTH_MULTIPLE = 1, and every next_multiple_of call degenerates to an identity.

The layout itself comes from one shared routine, pack_kgroup_panels in gemmkit/src/pack.rs. It is the single source of truth for the interleave index math.

Plain pack_panels stores a panel depth-major. At each depth step, it stores width contiguous leading elements: mr rows for the LHS, nr columns for the RHS.

pack_kgroup_panels instead groups the depth axis Q at a time, so one lane’s Q consecutive depth values become contiguous. Within a panel, group g, lane i, in-group position t lands at offset g*width*Q + i*Q + t.

That is exactly the shape one dot instruction reads. A 64-byte A register covers LANES rows times Q contiguous depth elements. A B group broadcasts Q contiguous depth values of one column, as a single 32-bit load.

2 more details of the shared packer matter. It takes a per-element transform, xform. That transform is the identity for bf16, and the +128 bias for VNNI’s A operand. The packer fills every pad position with xform(0). Pad positions are the leading positions past the block, and the depth positions past kc. This keeps the pad always consistent with the live elements.

The interleaved layout cannot be read in place. So every dot family sets FORCE_PACK_LHS and FORCE_PACK_RHS, overriding the driver’s cost-based pack decision. A dot kernel always pays the pack cost. This is precisely the cost the gates described below hedge against.

A byte-level oracle in pack.rs’s tests reimplements the layout naively. It checks that the real routine reproduces it bit-for-bit, across width tails, depth tails, and strided sources.

On the consuming side, KernelSimd::dot_accumulate is the seam the families call, instead of the widen-FMA loop. Its default body is unreachable!. Only a dot-capable token overrides it, and only a dot family ever calls it.

Avx512Vnni and Avx512Bf16 are distinct tokens from Avx512F, because #[target_feature] works per token. _mm512_dpbusd_epi32 needs an avx512vnni codegen context. Avx512F::vectorize only establishes avx512f, so it cannot provide that context.

The method receives the real, unpadded kc. It reads ceil(kc / Q) instruction groups from the depth-padded panels. Any signedness or bias correction is applied internally, so the accumulators hold the true sum_k(A*B) on return.

The fold lives on this dedicated seam, rather than on the generic accumulate_tile, for a documented reason. Fusing depth steps reshapes the accumulation rounding, and accumulate_tile’s contract forbids that.

i8 through vpdpbusd

vpdpbusd computes an unsigned-times-signed dot. It takes u8 from the first operand and i8 from the second. GEMM wants signed-times-signed instead.

The fix is algebraic, not per-element. The LHS pack offsets every byte by +128, into the unsigned domain. This transform is vnni_a_xform in gemmkit/src/kernel/int.rs. It uses the constant VNNI_A_BIAS = 128, defined once in gemmkit/src/simd.rs, so the pack and the correction can never drift apart.

sum_k((A+128)*B) = sum_k(A*B) + 128*sum_k(B). So the kernel recovers the true product by subtracting a per-column correction, 128 * sum_k(B[k][j]). Avx512Vnni::dot_accumulate computes those column sums with a small scalar pass over the signed packed B panel, before the vector loop. It then subtracts the splatted correction from every accumulator at the end.

The pads cooperate with this scheme. The A pad is xform(0) = 128, and the correction cancels its contribution exactly. The B pad is 0, so it contributes nothing to the product or to the column sums.

i32 accumulation wraps, and wrapping addition is associative modulo 2^32. So regrouping the sum into quads, plus the bias correction, equals the ascending-k widen sum bit-for-bit. IntGemmVnni and the widen IntGemm produce identical output on every input.

The ISA choice can therefore never change an i8 result. This is a stronger property than the reproducibility contract requires. That contract only promises reproducible results on a fixed machine and configuration, not bitwise agreement across kernel choices.

That freedom to swap kernels mid-flight is what the small-parallel fallback gate exploits. The VNNI pack is mandatory on both operands. On a small multi-threaded problem, that pack barrier can dominate the compute it is meant to save.

3 conditions trigger the fallback. The ISA selection must be automatic. The parallelism must be Rayon(n) with n != 1. And m*n*k must fall below GEMMKIT_I8_VNNI_MIN_PAR_MNK, whose default is 768^3. When all 3 hold, dispatch/int.rs hands the call to the in-place widen kernel instead.

Serial runs and large parallel runs keep VNNI. A forced GEMMKIT_REQUIRE_ISA=avx512vnni disables the gate entirely, because a forced pin must run exactly that kernel.

The prepacked-RHS path also bypasses the gate, for 2 reasons. A k-quad-interleaved buffer is only consumable by the VNNI family. The pack barrier the gate hedges against was already amortized once, at prepack time. VNNI’s RHS pack is otherwise mandatory on every call. So prepacking is a bigger win there than for any kernel that can read its operand in place.

bf16 through vdpbf16ps

Bf16DotGemm is the floating-point sibling. Its DEPTH_MULTIPLE is 2. Both operands are packed k-pair-interleaved, with each pair stored as one 32-bit __m512bh element. dot_accumulate issues one vdpbf16ps per accumulator, per pair-step.

Everything downstream of the accumulation is shared verbatim with MixedGemm<bf16>. This includes the alpha fold and the widen-read, narrow-store epilogue, through the common mixed_epilogue helper. The family keeps OUT_IS_ACC = false, so the whole contraction accumulates in f32 and rounds to bf16 exactly once.

The numeric story differs from VNNI in one essential way. vdpbf16ps’s fused 2-term dot rounds differently from 2 separate widen-FMAs. So the dot kernel’s result is only tolerance-equal to the widen path, not bitwise equal.

That is exactly what the engine’s consistency bar allows. Results must be reproducible under a fixed input, environment, and configuration. They need not be bitwise-identical across kernel choices. The dot kernel is fully deterministic. Serial, parallel, and prepacked runs all share the same kernel and pack layout, so they reproduce each other bit-for-bit.

There is no size gate on this path. Auto-selection prefers Bf16DotGemm whenever the CPU reports avx512bf16, because it is a structural win over the plain widen path. Unlike VNNI, there is no small-parallel fallback here.

3 special-path reroutes are the only exceptions. gemv, small_mn, and small-k shapes deliberately stay on MixedGemm<bf16>’s widen seam. A tiny or degenerate output folds nothing, and the dot pack’s depth padding is pure loss there. The i8 dispatch reroutes its own tiny shapes to the widen kernel, for the same reason.

The deep-K problem

OUT_IS_ACC = false buys single-rounding at a structural price. The driver runs kc = k, one depth panel spanning the entire contraction. This replaces the cache-model kc slices every homogeneous family gets (see Blocking and the Cache Model).

The RHS micropanel a microtile call reads is then nr * k * sizeof(N) bytes. Once that micropanel outgrows the L2 cache, every one of the m/mr microtile calls in a column strip streams it from L3 or DRAM instead. The even larger mr * k LHS micropanel streams from there too. The cliff this creates is sharp. Throughput stays near peak while the micropanel fits in L2, and falls once it does not.

The engage gate in gemmkit/src/dispatch/mixed.rs compares that micropanel size against a byte threshold:

#![allow(unused)]
fn main() {
let engage_deep_k = NR
    .checked_mul(t.k)
    .and_then(|x| x.checked_mul(core::mem::size_of::<N>()))
    .is_some_and(|bytes| bytes > crate::cache::deep_k_engage_bytes());
if engage_deep_k {
    run_deep_k_twin::<N, Fam::Twin, S, MR_REG, NR>(simd, &t, par, ws);
    return;
}
}

The threshold is the GEMMKIT_DEEP_KC_BYTES knob, taken verbatim when it is non-zero. The default, 0, derives the threshold as half the effective per-worker L2 capacity.

Half of L2, rather than the whole of it, is a deliberate choice. A gate set to the full L2 size would engage too late. It would fire well past the point where the micropanel no longer fits, and so it would miss the cliff. Half of L2 leaves room for the rest of the working set, and engages the twin while there is still time to avoid the cliff.

The checked_mul chain fails closed. A broadcast operand can pass validation with a logically absurd k. An overflowing size must fall through to the single panel instead. That panel’s own pack sizing then rejects the problem, rather than let a twin multi-slice that k forever.

The f32-output twin

Above the gate, dispatch does not run the narrow family at all. A small DeepKTwin trait maps each narrow family to its f32-output twin. MixedGemm<N> maps to MixedGemmF32<N>, and Bf16DotGemm maps to Bf16DotGemmF32. The only type change in each twin is Out = f32 = Acc.

That one change flips OUT_IS_ACC back to its default of true. The driver’s ordinary multi-slice K blocking then applies unchanged, and every slice’s panels are L2-resident again. That is the entire point of the twin.

The pack layout and the accumulation loop are the narrow family’s, verbatim. MixedGemmF32 reuses pack_panels and the shared widen-FMA helper. Bf16DotGemmF32 reuses pack_kgroup_panels and dot_accumulate. The accumulation helpers only touch the input side of the KernelSimd seam. So the accumulator they produce is byte-for-byte what the narrow family would compute.

The twin runs with alpha = 1 and beta = 0, into an m x n column-major f32 scratch buffer. That buffer is drawn from a dedicated Workspace. Deep-K is by definition a large-k regime, so one m*n f32 allocation is negligible. Keeping it separate also leaves the pooled packing workspace free for the twin driver to use.

Afterward, one vectorized sweep computes narrow(alpha*scratch + beta*widen(C)). This replicates mixed_epilogue’s arithmetic operation for operation, including the same store_out narrowing.

What makes the route more than an approximation is how the slices chain together. The twin’s microkernel seeds its accumulator registers from the scratch buffer, through a third KernelSimd<N, N, f32, f32> seam, twin_seed in gemmkit/src/kernel/mixed.rs. On an accumulate slice, it loads the running partial into the registers and continues the ascending-k chain. It never sums a slice from zero and adds the result afterward.

A store and reload of an f32 is exact. So the multi-slice sum is exactly the single-panel sum, merely split at slice boundaries. For beta in {0, 1}, the deep-K result is byte-for-byte the single-panel result.

For a general beta, the result holds only to tolerance instead. The reason is mundane. The single panel fuses beta*C + AB in one FMA on full tiles, but combines the same terms unfused on edge tiles. No single sweep can match both cases at once.

Serial and parallel runs remain bit-identical in every case. The twin driver’s blocking does not depend on the thread count, and the final sweep is elementwise.

The dot twin has one extra alignment rule. The driver rounds the blocking kc up to DEPTH_MULTIPLE, so an interior slice boundary never splits a k-pair. A split pair would zero-pad mid-contraction and regroup the fused dot incorrectly. With this rule, only the final short tail is ever padded, exactly as in the single-panel case.

3 routes deliberately keep the single panel. They are shallow k below the gate, the fused-epilogue entries, and the prepacked-RHS path. On the prepacked path, a DEPTH_MULTIPLE > 1 buffer requires the whole contraction to stay one depth slice. The driver enforces this with a hard assert, because violating it would misaddress micropanels silently.

The parity claims are tested directly. gemmkit/tests/deep_k_narrow.rs toggles GEMMKIT_DEEP_KC_BYTES. A value of 1 forces the twin at any k, and usize::MAX forces the single panel. The test checks byte-for-byte equality for beta in {0, 1}, and tolerance for general beta, on whichever ISA the host selects. Tuning Knobs documents this knob along with the rest.

The Complex Split Kernel

Complex GEMM is the one homogeneous type family that does not ride FloatGemm. On paper it could. Complex<f32> and Complex<f64> accumulate in themselves, so Lhs = Rhs = Acc = Out. That is exactly the shape the float family handles.

The problem is the memory layout. num_complex stores a complex number as an adjacent (re, im) pair. So a SIMD register loaded from a complex slice holds re, im, re, im, .... One complex multiply needs cross-lane combinations: re*re - im*im for the real part, and re*im + im*re for the imaginary part.

On interleaved lanes, those combinations force shuffle and fmaddsub-style instructions inside the innermost loop. That cost repeats once per depth step, O(mnk) times in total. gemmkit instead rewrites the layout once, at pack time, and keeps the hot loop as pure real FMAs, with no in-loop shuffles at all.

This page covers 5 things. First, the split design itself. Second, how conjugation falls out of it for free. Third, the seam the kernel runs through. Fourth, the register-budget arithmetic behind the tile shapes. Fifth, what the numerics guarantee.

The split layout

The family is ComplexGemm<T, CONJ_A, CONJ_B>, in gemmkit/src/kernel/complex.rs. Its packers are where the design lives.

pack_planar lays each micropanel down in structure-of-arrays form. For every depth step, the panel stores width real parts, immediately followed by width imaginary parts. width is mr for the LHS and nr for the RHS, with strides swapped exactly as in the shared pack_panels.

The kernel then loads a register of reals and a register of imaginary parts with plain contiguous loads. The de-interleave cost moves out of the kc inner loop and into the pack step. The amortized cost becomes O(MK + KN), instead of O(MNK).

The kernel can only consume that planar layout. So both operands are always packed. The family sets FORCE_PACK_LHS = FORCE_PACK_RHS = true, overriding the driver’s cost-based decision to read an operand in place instead.

pack_planar mirrors pack_panels’ 2 write paths. One is a straight walk, used when the leading dimension is contiguous. The other is a cache-blocked transpose, used for a strided source, so a row-major operand packs without a cache miss per element. The 2 paths write byte-identical panels. Only the write order differs. See Packing and Workspaces for the shared framework.

Conjugation as a pack-time sign flip

Conjugation only negates the imaginary part. The pack already writes the imaginary plane separately. So conj(A)*B and A*conj(B) cost nothing in the hot loop.

Setting the CONJ_A or CONJ_B const generic makes the packer negate the imaginary plane as it copies. This is a true negation, so +0.0 maps to -0.0, matching num_complex’s .conj(). The same real-FMA loop then runs unchanged, with no per-element conjugation branch anywhere.

This is also the second reason for the force-pack flags. When packing does more than a plain copy, the transform must always run.

The runtime-to-compile-time bridge lives in gemmkit/src/dispatch/complex.rs. The public entry, gemm_cplx, takes conj_a and conj_b as plain bools. run_complex matches on that pair once, and dispatches to the matching one of 4 ComplexGemm monomorphizations. The branch happens once per call, never inside the loop.

One subtlety rides along with this. The orientation swap that normalizes a row-major-ish C computes C^T = B^T * A^T instead. Since (conj(A)*B)^T = B^T * conj(A)^T, the swap must also swap the conj flags. It does this before the match runs.

Output conjugation, conjC, is not implemented. On the degenerate path, where k == 0 or alpha == 0, the flags are simply irrelevant. With no A*B term, there is nothing to conjugate.

The hot loop: 4 real FMAs per complex MAC

There is one more layering problem to solve before the loop can run. The family is homogeneous, so the driver’s bound is KernelSimd<T, T, T, T> with T complex. That bound only yields SimdOps<Complex<..>>, not the real-typed operations the split kernel actually needs.

The bridge is the SimdOps::cplx_microkernel seam. The family’s microkernel forwards to it. Each ISA token’s override, generated by the impl_complex_simd! glue macro in gemmkit/src/simd/complex.rs, forwards in turn to one shared, ISA-generic function, soa_microkernel, written over S: SimdOps<C::Real>.

The accumulator stays Complex-typed at the family seam. So complex alpha and beta thread through the driver unchanged. Inside the seam, though, the accumulators are 2 banks of real registers.

The thin SimdOps<Complex<..>> glue exists only so the driver can read LANES, and so the homogeneous blanket impl applies. Its element operations are all unreachable!, because complex GEMM never calls them. LANES is set to the real lane count. So one real lane spans one complex row, and the driver’s mr = MR_REG * LANES counts complex rows.

The loop itself, from gemmkit/src/simd/complex.rs:

#![allow(unused)]
fn main() {
for p in 0..kc {
    let are_p = a_re.add(p * 2 * mr); // re plane of this depth step
    let aim_p = are_p.add(mr); // im plane (offset by `mr`)
    let ar: [<S as SimdOps<C::Real>>::Reg; MR_REG] =
        core::array::from_fn(|i| simd.loadu(are_p.add(i * lanes)));
    let ai: [<S as SimdOps<C::Real>>::Reg; MR_REG] =
        core::array::from_fn(|i| simd.loadu(aim_p.add(i * lanes)));
    let bre_p = b_re.add(p * 2 * NR);
    let bim_p = bre_p.add(NR);
    for j in 0..NR {
        let br = simd.splat(*bre_p.add(j));
        let bi = simd.splat(*bim_p.add(j));
        for i in 0..MR_REG {
            acc_re[j][i] = simd.mul_add(ar[i], br, acc_re[j][i]); // += ar*br
            acc_re[j][i] = simd.fnma(ai[i], bi, acc_re[j][i]); //    -= ai*bi
            acc_im[j][i] = simd.mul_add(ar[i], bi, acc_im[j][i]); // += ar*bi
            acc_im[j][i] = simd.mul_add(ai[i], br, acc_im[j][i]); // += ai*br
        }
    }
}
}

One complex multiply-accumulate is 4 fused real steps, into 2 running banks. acc_re gets a mul_add and an fnma, the fused negate-multiply-add, vfnmadd on x86. acc_im gets 2 mul_adds.

Every operation is a plain lane-parallel FMA, on contiguous loads and scalar broadcasts. Nothing crosses lanes until the epilogue. The fixed per-p order is deliberate. It is what makes full tiles and edge tiles of the same matrix round identically.

After the loop, the banks drain to planar scratch. A scalar epilogue folds complex alpha, skipping the complex multiply when alpha == 1. It combines beta*C case by case, and re-interleaves on store. This is an amortized O(MN) pass that handles full, edge, and strided output tiles uniformly.

The scalar de-interleave in the pack, and the scalar re-interleave in the epilogue, are both a deliberate choice, not an oversight. The inner loop dominates the total cost of the kernel. So the generic scalar path stays the floor for both steps, on every ISA.

Register pressure and the NR choice

The split design doubles the accumulator count. An MR_REG x NR complex tile needs 2*MR_REG*NR accumulator registers, a real bank and an imaginary bank. It also needs 2*MR_REG A-plane registers, and 2 B splats per column step. That budget is documented tile by tile in gemmkit/src/dispatch/complex.rs. It makes the complex tiles smaller than the float ones, and it is where the tile shapes on this page come from.

On FMA, with 16 YMM registers, c32 runs MR_REG = 1, which is 8 complex rows at 8 real lanes, with NR = 5. That is 10 accumulators, plus 2 A registers, plus 2 B splats, for 14 of 16 registers. The 2 spare registers matter. A full 16-of-16 tile with NR = 6 instead spills accumulators to the stack. So NR was shrunk to 5, and the code comment records why.

AVX-512’s 32 ZMM registers relax the pressure. c32 runs MR_REG = 2, NR = 6, using 24 + 4 + 2 = 30 of 32 registers. NEON, with 32 vector registers, runs MR_REG = 2, NR = 5, for 26 of 32, leaving room for in-flight load temporaries. wasm simd128 runs MR_REG = 1, NR = 4, for 12 live v128 registers.

Each c64 variant keeps its c32 sibling’s MR_REG and NR, with halved lanes. The budget arithmetic does not depend on the lane count.

Accuracy and reproducibility

Complex has no special paths. run_complex routes every shape to driver::run, with no gemv, small_mn, or small-k arms. The Special Paths machinery is real-float and integer only.

This makes the numeric contract easy to state. gemmkit/tests/correctness/complex.rs tests each piece of it directly.

Determinism and thread independence hold bitwise. Blocking is thread-count independent, so serial and parallel runs of the same problem produce bit-identical output. The test asserts equality of the raw re/im bit patterns across thread counts.

Within one run, the fixed per-step order of the 4 FMAs means full tiles and edge tiles of the same matrix round identically. So results do not depend on where a tile boundary happened to fall.

Conjugation adds no rounding at all. It is a sign flip on exact values. A dedicated test on small-integer inputs, where every product and sum is exactly representable, checks all 4 conjugation combinations against a naive reference. That test uses exact equality, not tolerance.

Against an external oracle, the bar is necessarily looser. The correctness suite compares gemm_cplx, including every conjugation combination, against the gemm crate, under the suite’s usual L2-style tolerance. It separately exercises a negative-row-stride view, plus conjugation, against a row-reversed reference. This looser bar exists because a blocked SoA contraction legitimately rounds differently from another engine’s ordering.

That is the reproducibility contract, applied to complex numbers. Under a fixed input, environment, and configuration, the result is the same bits. Across different engines or different orderings, results agree only to tolerance.

The fused-bias entry is the exception that gets a bitwise guarantee. gemm_cplx_fused supports a per-row or per-column complex bias. It deliberately supports no activation, since ReLU-style activations are undefined on an unordered field.

Its epilogue never touches the kernel’s arithmetic. The SoA kernel stores exactly the bits plain gemm_cplx would store. A tile-local post-pass then maps them in place, on the final depth panel only. Complex is OUT_IS_ACC = true, so intermediate panels must keep their raw partials.

The result is bit-identical to gemm_cplx followed by the same element-wise bias add, for every shape and every conjugation combination. The Identity instantiation const-folds the post-pass away entirely. So the non-fused path pays nothing for the hook’s existence. Epilogue Fusion covers the general mechanism.

Blocking and the Cache Model

The microkernel at the bottom of the driver multiplies one MR x NR tile of C. It holds the tile in registers and streams a micropanel of A and a micropanel of B along the depth axis. This loop reaches machine peak only when both streams come from nearby cache. A GEMM touches far more data than any cache holds. It reuses every element of A across n output columns, and every element of B across m output rows.

Blocking is how the driver arranges this reuse. It partitions the problem so each operand block enters a specific cache level once, then gets read many times before anything evicts it. KC slices the depth axis so the 2 micropanels a tile multiplication reads stay resident in L1 for the whole tile. MC sizes the packed A macro-panel so it stays in L2 while the driver sweeps it across every column tile of the current column block. NC sizes the packed B macro-panel so it stays in L3 while every row block sweeps over it. Panel residency across the loop nest is the entire point. Without it, the same bytes would stream from DRAM m, n, or k times over.

Many libraries hard-code (MC, KC, NC) per microarchitecture. gemmkit instead computes them analytically for each call, in CacheTopology::blocking (gemmkit/src/cache.rs, layer L3). The function follows the BLIS model and derives the sizes from cache geometry detected at runtime. Its inputs are the microtile (mr, nr), the byte size of one packed input element, and the problem shape (m, n, k). Its output is the Blocking { mc, kc, nc } triple that the driver’s loop nest iterates by. Life of a GEMM Call shows where each value lands in the nest. This page explains how the model derives the values, and why it takes this shape.

3 constraints, 3 block sizes

KC: both micropanels in L1 without self-eviction

Every microtile call walks kc depth steps. Each step reads mr packed A elements and nr packed B elements. So an mr x kc micropanel and an nr x kc micropanel must both fit in L1d for the whole tile. The subtle requirement is without self-eviction. A cache is not a byte pool. It holds sets of ways. A panel that maps too many of its own lines onto the same sets evicts itself before its total size even reaches the cache size.

The model therefore works in lines and sets, not bytes. It computes how many L1 lines one depth step of each micropanel claims. Then it picks the largest kc whose combined footprint stays within the L1 associativity. It raises that result to the GEMMKIT_KC_MIN floor (default 512), so a small L1 never starves the microkernel’s depth walk, and clamps it to k. A final rebalance splits k into ceil(k / kc) panels of near-equal size, so the last depth slice is never a sliver.

MC: the A macro-panel in L2, minus what B needs

Within one row block, the driver reuses the packed mc x kc A panel across every column tile of the current column block. So the panel should fill L2. It cannot fill all of L2, though. The nr x kc B micropanel of the same depth slice also streams through L2 on every tile call. The model counts how many L2 ways that micropanel occupies. It reserves those ways plus one spare way, then hands the rest of the capacity to A.

It divides that remaining capacity by kc to get mc, and rounds the result down to a multiple of mr. It then rebalances the result so the row blocks come out even. A BLIS-style hard cap of GEMMKIT_MC_REG_PANELS * MR rows (default 8 microtile rows) clamps the result last. This cap is a calibration point, not a bound strictly derived from the L2 term. In practice the cap binds before the L2 capacity term does. So MC ends up as a small multiple of MR, and the L2 capacity term mostly serves as headroom.

NC: the B macro-panel in L3, or a panel cap without one

When an L3 is present, the model reserves one way for the A traffic passing through. It budgets the rest for the packed kc x nc B macro-panel. It divides that capacity by kc to get nc, rounds the result down to a multiple of nr, and rebalances it across n.

Some machines report no L3 at all. On Apple Silicon, for example, a cluster-shared L2 tops the hierarchy. There, the model runs full-N up to a panel-count cap instead. nc is GEMMKIT_NC_NO_L3_PANELS * nr (default 512 panels, that is, 2048 columns at nr = 4), capped by n. With no L3 to keep B resident, B streams from DRAM regardless. The cap only bounds the shared packed-B buffer. It does not model residency.

Sized in packed elements, not accumulator elements

The sizeof argument is the size of one packed input element. The driver passes size_of::<Fam::Lhs>(), not the accumulator size, because the model stores the panels it budgets in packed Lhs/Rhs units. For f32 and f64 the 2 sizes coincide, so nothing changes.

For narrow types the distinction matters more. i8 packs 1 byte per element against a 4-byte i32 accumulator. f16 and bf16 pack 2 bytes against a 4-byte f32 accumulator. Sizing by the accumulator would cut their kc and nc to a quarter or half of what the caches actually fit. Narrow types instead get proportionally deeper blocks, which is why they can outperform f32 on the same hardware. The prepack entries reuse the same model with a sentinel row count, so a prepacked operand’s geometry stays independent of the eventual m.

Prefetching the output tile past the LLC

The 3 block sizes keep the A and B panels resident, but they say nothing about C. Every microtile call reads, modifies, and writes its mr x nr output tile. Once a call’s working set (the A, B, and C bytes together) outgrows the per-core-reachable LLC, that output tile no longer lives in cache. Its store then reaches into DRAM.

The driver answers with a software prefetch, decided once per call. It compares the working set against cache::prefetch_ws_bytes (the GEMMKIT_PREFETCH_MIN_BYTES gate, where 0 means auto: the per-core-reachable LLC, L3 where present, otherwise L2). Once the gate clears, the driver issues a T0 prefetch of each output microtile just ahead of that tile’s microkernel call. This pulls the lines the store will touch into L1 while the microkernel still computes.

The hint walks whole 64-byte cache lines along the tile’s unit-stride dimension. A tile strided in both dimensions has no contiguous lines, so the driver skips it. The driver emits the hint only on x86_64 (prefetcht0, baseline SSE, needs no feature gate). On any other target it lowers to nothing, so aarch64 and wasm stay unaffected. The prefetch moves cache lines only, never arithmetic, so results stay bit-identical whether the gate is on, off, or forced. Below the gate, where the tiles stay cache-resident, the prefetch path adds no extra cost.

The tiny-matrix shortcut

When both m and n are at or below GEMMKIT_TINY_BLOCK_DIM (default 64), the driver skips the full model. It sets kc to k, clamped to the GEMMKIT_KC ceiling (default 2048, or 16384 on aarch64). It sets mc to whatever row count keeps the panel in L2 at that depth, capped by m itself. It sets nc to n rounded up to nr.

A problem whose whole working set fits in L2 gains nothing from 3 levels of residency analysis. The shortcut spends the saved arithmetic where it matters: on the fixed per-call overhead that dominates small products.

The ceiling counts 4-byte elements, and a narrow element divides it. So f16 gets 2 times the depth of f32, and int8 gets 4 times. What stays fixed is the byte budget, which is what the hardware limit is about. One number then calibrates every element family, and a tuner sweep run on f32 transfers to the rest.

The ceiling sets the depth-slice count, and the slice count drives 2 costs that pull in opposite directions. Each extra slice re-reads and re-writes C, re-enters the driver, and forks the workers once more. That argues for a deep ceiling.

A deeper slice also grows the packed A and B panels. Those panels must stay in a private L2, so that argues for a shallow ceiling. On x86 the panels reach about 1.1 MiB at the default, which is one Zen5 L2. The parallel cost is the larger of the 2, so the default sits at the residency limit rather than below it.

A wide element is the case where those 2 costs disagree about the divisor. A byte budget is not a slice budget. At a fixed budget, a 16-byte element takes 4 times the slices of a 4-byte element, so it pays the per-slice cost 4 times as often. Residency wants the division, and the slice count does not.

Which cost wins is a property of the machine, so the divisor carries an arch-split cap. On x86 the private L2 is 1 MiB, residency binds first, and the divisor applies at every element size. On aarch64 the effective L2 is 4 times larger, and unified memory feeds an overflowing panel well. Residency barely binds there, so a wide element keeps the whole ceiling. Measured on an M4 Max, dividing it cost c64 21 to 51 percent on the parallel path.

Detection: a fallback chain that cannot fail

The model is only as good as the geometry it receives, and there is no portable way to ask for cache geometry. gemmkit instead runs a best-effort chain, where #[cfg] only ever picks the sniffing method, never the values. A #[cfg(target_arch)] check cannot tell an Intel part apart from an AMD part, and a VM or container can mask CPUID or hide /sys. So every backend returns an Option, and the chain bottoms out in a constant that cannot fail.

#![allow(unused)]
fn main() {
// gemmkit/src/cache.rs
#[cfg(feature = "std")]
fn detect() -> CacheTopology {
    // try the CPUID backend
    #[cfg(all(any(target_arch = "x86", target_arch = "x86_64"), not(miri)))]
    if let Some(t) = cpuid::detect().filter(plausible) {
        return t;
    }
    // try the sysfs backend
    #[cfg(all(target_os = "linux", not(miri)))]
    if let Some(t) = sysfs::detect().filter(plausible) {
        return t;
    }
    // try the sysctl backend
    #[cfg(all(target_os = "macos", not(miri)))]
    if let Some(t) = sysctl::detect().filter(plausible) {
        return t;
    }
    ZEN5_FALLBACK
}
}

The backends run in this order.

CPUID (cache/cpuid.rs) is a single instruction, so it works regardless of OS, in containers, and in most VMs. CPUID reads both vendors through the per-cache topology leaf (Intel 04h, AMD 0x8000_001D), which describes each cache as reachable from the executing core. On a multi-die part such as a 2-CCD Ryzen, the L3 figure is the one complex a core can reach (32 MiB on the 9950X). It is not the package total. That per-core-reachable figure is the semantic every consumer of the value wants. AMD parts or hypervisors without that leaf fall back to the legacy L1 (0x8000_0005) and L2/L3 (0x8000_0006) leaves. There, L3 size arrives in units of 512 KiB as a package total, and the legacy fields cannot even encode 16-way associativity.

Linux sysfs (cache/sysfs.rs) parses /sys/devices/system/cpu/cpu0/cache/index*/ with plain std::fs. It serves as a fallback on x86 Linux, for a hypervisor that masks CPUID. It is also the primary source on aarch64 Linux, which has no CPUID instruction.

macOS sysctl (cache/sysctl.rs) reads sysctlbyname keys through a 2-line extern "C" block, with no libc dependency. It prefers the Apple Silicon per-performance-level keys (hw.perflevel0.*, the P-cores), with the flat Intel-Mac keys as a fallback. sysctl does not expose associativity, so the backend assumes conservative typical values. This is safe, because the model clamps associativity with .max(2) and needs it only approximately.

The bottom of the chain is ZEN5_FALLBACK, a static default calibrated on the Ryzen 9950X dev machine. L1d is 48 KiB and 12-way. L2 is 1 MiB, 16-way, and private. L3 is 32 MiB and 16-way.

2 guards make the chain robust rather than merely ordered. plausible rejects half-populated reads. Any level smaller than 4 KiB, a line under 16 bytes, or zero associativity fails the whole backend. A masked leaf therefore cannot poison blocking with zeros. Detection also runs at most once per process. Machine::current() memoizes the topology behind a OnceLock. It also memoizes the OS page size (getpagesize, validated as a power of two between 4 KiB and 2 MiB). That page size drives the LHS-packing stride gate described in Packing and Workspaces. A no_std build skips detection entirely and uses the Zen5 fallback with a 4 KiB page.

shared_by: contention for what the driver puts there

Each Level carries bytes, assoc, line, and one derived field, shared_by, which divides the level’s capacity into the effective_bytes that the model actually budgets. Storing the hardware core-sharing count there would seem natural, but it would be wrong. shared_by instead models per-worker contention for the data the driver actually places at that level. The driver’s placement is: per-worker A/B micropanels in L1d, each worker’s private A macro-panel in L2, and one shared B macro-panel in L3.

That placement fixes the values. L1d is per-core, so its whole capacity serves one worker’s micropanels, and shared_by = 1. L3 is shared by every core in hardware. The data the driver keeps there, though, is a single panel. All workers read that panel in common: the same bytes, not per-worker copies. So the whole level belongs to that one panel, and shared_by is again 1. Dividing by the raw core count would shrink the budget many times over and crater NC for no reason.

Only L2 holds genuinely private per-worker data, so only L2 uses the physical-core L2-sharing degree. That degree is 1 on parts with a private L2, such as mainstream x86 and Neoverse. It is the cluster size on parts where a core cluster shares one L2, such as Apple Silicon. There, several workers’ private A panels really do contend for the same ways. Each backend must derive this value rather than copy a raw count. sysfs divides the raw L2 shared_cpu_list count by the SMT degree read from L1d’s sharing list, so hyperthread siblings are not double-counted. sysctl reads hw.perflevel0.cpusperl2. The CPUID backend hard-sets 1, because x86 L2s are private to each physical core.

On x86 and Graviton the whole mechanism reduces to all-ones. It exists for the cluster-L2 parts. There, it decides whether the model blocks for the L2 a worker actually gets, or for one it must share across a whole cluster.

What the thread count moves, and what it cannot

blocking takes no thread-count parameter. For KC and NC that omission carries real weight. Both depend only on the machine and the problem. So a serial run and a wide parallel run derive the same KC and the same NC. That gives every run the same depth slices, and the same fixed-order depth chain for every output element.

MC is the one blocking dimension the driver does adjust for parallelism. A wide worker count can leave the flat job list too shallow, with fewer than a handful of chunks per worker. The run’s tail then degenerates into idle workers waiting on whoever drew the last chunks. When that happens, the driver shrinks MC to cut more row blocks and deepen the list. Parallel Execution details this parallel job-depth floor. So the panel boundaries and the flat job list are not strictly worker-count independent anymore.

Bit-identity survives the MC shrink regardless, because the shrink itself carries no numerics. MC always stays an MR multiple. So the set of microtiles it produces (every MR-aligned row offset plus the single m-tail tile) stays identical under any split. A wider worker count only regroups the same tiles into more, smaller row blocks. KC, the only blocking dimension that shapes a tile’s accumulation order, never moves with the thread count. So under a fixed configuration, changing only the worker count leaves every output element’s accumulation order unchanged. That is the mechanism behind gemmkit’s reproducibility contract. Parallel Execution assembles the full contract and states its exact scope.

Parallelism otherwise influences packing decisions only. The LHS pack gate reads per-worker column reuse, and the shared-A pre-pass engages only on large parallel problems. Those decisions choose where the driver stages packed bytes and who writes them. They never change what values the kernel computes. Parallel Execution covers how the job list splits, and how the contract holds end to end. Tuning Knobs catalogs every GEMMKIT_* threshold named on this page, and every other one.

Packing and Workspaces

The microkernel wants its inputs in exactly one shape. For each depth step, it wants mr elements of A contiguous in memory, and nr elements of B contiguous in memory. It wants them panel after panel, with nothing between them. User matrices almost never look like that. They have arbitrary row and column strides, tails that do not divide the microtile, and depth walks that can stride across memory pages.

Packing is the copy that closes this gap. It rearranges each macro-block into micropanel-major layout once, so the innermost loop reads pure unit-stride streams, full mr/nr vectors every time, from 64-byte-aligned scratch. The copy costs O(mc*kc), against the O(mc*kc*nc) compute that reuses it. That is why it amortizes, and why the driver still skips it when reuse is too low to pay for it.

One routine, both operands

The mechanical copy lives in a single routine, pack_panels (gemmkit/src/pack.rs, layer L1). The LHS and RHS layouts are the same layout, viewed from different sides. An LHS macro-block packs into panels mr rows tall, stored column by column. Panel 0 holds rows 0..mr, with each depth step’s mr elements contiguous. Panel 1 holds rows mr..2*mr, and so on. An RHS macro-block packs into panels nr columns wide, stored row by row.

Both layouts are width contiguous leading elements per depth step. The only difference is which matrix axis plays the leading role. So the 2 KernelFamily hooks call the same routine, with the strides swapped:

#![allow(unused)]
fn main() {
// gemmkit/src/kernel/float.rs
#[inline]
unsafe fn pack_rhs(
    dst: *mut T,
    src: *const T,
    rs: isize,
    cs: isize,
    kc: usize,
    nc: usize,
    nr: usize,
) {
    // RHS panels are `nr` columns wide, stored row-by-row: the "leading"
    // direction is columns (stride `cs`) and the "depth" is rows (stride
    // `rs`), the transpose of the LHS case, handled by swapping strides
    unsafe {
        pack_panels(
            dst, src, /*lead*/ cs, /*depth*/ rs, /*n_lead*/ nc, kc, nr,
        )
    }
}
}

pack_lhs is the mirror image: lead = rs, depth = cs, width = mr. When the block does not divide evenly, the routine zero-fills the tail panel’s dead lanes. The kernel then always reads full mr/nr vectors, and edge tiles need no masking in the multiply itself.

Inside the routine, 2 paths write byte-identical output. The 1st case is a contiguous leading dimension (lead == 1), for a column-major A or a row-major B. There, each depth step’s live elements already sit adjacent in the source. The panel is then a sequence of straight copy_nonoverlapping calls, plus a tail zero-fill.

The 2nd case is a strided leading dimension. There, a naive gather would take a cache miss per element (width strided loads per depth step). Instead the routine runs a cache-blocked transpose. It walks the source along its contiguous dimension in strips of GEMMKIT_PACK_TRANSPOSE_TILE depth steps (default 16), and scatters each strip into the panel. This produces the same packed bytes as a pure reordered copy, but far more cheaply for a strided source. That is what makes row-major-A layouts cost little more than column-major ones.

The dot-product families (i8 VNNI, bf16 vdpbf16ps) have a sibling routine, pack_kgroup_panels. It additionally interleaves DEPTH_MULTIPLE consecutive depth steps per lane, so one dot instruction can consume a whole group. That layout belongs to Dot Kernels and the Deep-K Twin.

Whether to pack at all is the driver’s call, and the 2 operands are asymmetric.

The microkernel reads A as mr-wide vectors. So the driver must pack A whenever its rows are not unit-stride, or the row panel is partial. Beyond that, the driver packs A when each worker’s column reuse clears GEMMKIT_LHS_PACK_THRESHOLD (default 256 columns on aarch64, 1024 elsewhere).

The driver also packs a column-major A under one more condition. Its depth walk must be page-scale in stride, wide in span, and reused by enough column tiles to be worth the cost. All 3 conditions must hold together:

1. The per-step stride reaches half a memory page (GEMMKIT_LHS_PACK_STRIDE, auto-derived from the page size memoized in Machine).
2. The whole depth-slice walk (csa * sizeof(Lhs) * kc) reaches GEMMKIT_LHS_PACK_SPAN bytes (auto: 4 MiB).
3. At least GEMMKIT_LHS_PACK_REUSE nr-wide column tiles reuse each packed panel (min(n, nc) / nr, rounded up, default 128 on x86, 4 on aarch64).

Each gate rules out a different case where packing would not pay for itself. A page-scale stride over a span that still fits in cache only re-walks lines that are already warm. Reading A in place then costs less than paying for a pack. The span condition keeps A in place until the walk is wide enough to thrash the TLB, regardless of how much reuse follows.

The reuse floor prices the opposite failure mode. Take a tall, skinny shape, where m is much greater than n. It reaches a large span from very few column tiles. Packing it would then amortize an expensive copy over too little reuse to be worth it.

The reuse floor differs by architecture, because the pack-versus-read-in-place trade differs by architecture. On x86 a pack costs more relative to an in-place read. So the driver waits for more reuse before it pays for one, and the default floor is 128 tiles. On aarch64 a pack costs less relative to an in-place strided read. So the driver packs sooner, and the default floor is 4 tiles.

B, by contrast, is only ever read by broadcasting single elements, so any layout works unpacked. The driver packs B purely for reuse: once per depth slice, when m clears GEMMKIT_RHS_PACK_THRESHOLD (default 2048) and enough row blocks will re-read it. Who performs these packs, and the barriers between packing and compute, are scheduling questions. Parallel Execution covers them.

Prepacked operands

A per-call pack is wasted work when the same operand appears in call after call. This is the inference pattern: fixed weights multiplied against a stream of activations. The prepack entries in gemmkit/src/api/packed.rs pack a whole operand once, up front. prepack_rhs walks any-layout B through its strides and returns a PackedRhs<T>. gemm_packed_b then multiplies against it, skipping the per-call RHS pack entirely. Prepacked Operands covers the usage side of this API. Architecturally, 3 properties matter.

First, the buffer records the blocking geometry it was built for: nr, kc, and nc. The consuming call reads that geometry back verbatim. The driver substitutes the recorded kc and nc for its own model output. Only mc still derives from the real m. So panel addresses always match what the buffer packed, even if a tuning knob changed between the pack call and the consuming call. The geometry itself resolves through the same blocking model as a plain call, with a sentinel row count of tiny_block_dim() + 1. This keeps it off the tiny-matrix branch, so it stays independent of the eventual m.

Second, the layout has one source of truth. prepack_rhs fills the buffer through driver::pack_rhs_full. This lays panels down in exactly the order the driver’s own per-slice pack writes them. The order is jc blocks outermost, then depth slices, then the nr-wide panels of each slice. The prepacked bytes therefore equal the per-call packed bytes. So a prepacked GEMM reproduces a plain gemm under the same configuration. The documented exceptions are tiny products (m and n at or below tiny_block_dim) and gemv-shaped products. Plain gemm reroutes those to special paths, so they may differ in the last ULP.

Third, the buffer stays read-only during the GEMM, so every worker shares it with no synchronization. Unlike the per-call B pack, it needs no barrier.

PackedLhs costs almost no extra code, because of the engine’s A/B symmetry. An m x k LHS is the RHS of the transposed product C^T = B^T*A^T. So prepack_lhs delegates to prepack_rhs_unchecked with the strides swapped, and gemm_packed_a consumes it through the transposed problem. This symmetry also explains the orientation asserts. A prepacked B requires a column-major-ish C (|csc| >= |rsc|), and a prepacked A requires a row-major-ish one. The other orientation would make dispatch swap the operand roles, and the baked-in layout could not serve that swap.

The int8 feature adds a heterogeneous twin, prepack_rhs_i8 and gemm_i8_packed_b, with 3 deliberate differences.

First, prepack_rhs_i8 pins its layout to whichever integer kernel the process’s memoized dispatch selected: the VNNI k-quad-interleaved layout, or the widen kernel’s plain panels. The consuming entry always runs that same family, so it can never misread the buffer.

Second, it rounds the buffer depth up to the dot kernel’s DEPTH_MULTIPLE = 4 and packs the whole contraction as one depth slice. This satisfies the driver’s single-slice guard for depth-padded families.

Third, it deliberately bypasses the dynamic small-parallel widen fallback that plain gemm_i8 applies below GEMMKIT_I8_VNNI_MIN_PAR_MNK. A vpdpbusd buffer is quad-interleaved, so the widen kernel simply cannot consume it. Because integer accumulation is exact, the result is bit-identical to plain gemm_i8 either way.

Prepacking matters most on exactly this path. The VNNI RHS pack is otherwise mandatory on every call, so at small m the per-call O(k*n) pack cost dominates the O(m*k*n) compute.

The workspace

All of this packing needs scratch memory. Workspace (gemmkit/src/workspace.rs) is its allocator: a growable buffer, 64-byte aligned (enough for AVX-512 stores), that grows to the next power of 2 and never shrinks. Per call, Workspace::regions carves it into a_regions equal LHS regions plus one shared RHS region, with each region rounded up to the alignment.

The LHS region count is the worker count on the per-worker pack path, or the row-block count on the shared-A path. The carving works the same way either way. When neither operand packs, the driver skips the reservation entirely, so an all-in-place workload never grows the pool.

Fail closed at the byte product

The sizing arithmetic is where a memory-safety subtlety hides. gemmkit accepts broadcast (zero-stride) views. These pass bounds validation with a tiny backing slice, while presenting logical dimensions up to isize::MAX. So the products that size the pack buffer can genuinely overflow usize. A wrapped (too-small) size would then under-allocate a buffer that the pack writes past.

The driver guards its element-count products with checked_mul, but element counts alone are not enough. Take k = 2^56 on the mixed-precision path, where kc == k. An LHS region of mc * kc elements, say 32 * 2^56 = 2^61, fits usize comfortably and sails through every element-level check. Multiply that count by the element size and round up to the 64-byte alignment, though, and the value wraps. The overflow only appears at the element-to-byte conversion. That is where the guard must sit: at the chokepoint every region size funnels through:

#![allow(unused)]
fn main() {
// gemmkit/src/workspace.rs
fn region_bytes(elems: usize, esize: usize) -> usize {
    elems
        .checked_mul(esize)
        .and_then(|b| b.checked_next_multiple_of(ALIGN))
        .unwrap_or_else(|| workspace_too_large())
}
}

Workspace checks every step: the byte product, the alignment round-up, the region sum, and the final A + B total. Any overflow panics with the same “too large” contract as the driver’s own sizing. This is fail closed. The code rejects an absurd problem loudly, instead of corrupting memory. The driver runs its element-count guards unconditionally for the same reason, even on routes that end up packing nothing. Skipping them would also skip the abort, and send the absurd k into the in-place loops to spin for an effectively unbounded time.

The pool, _with, and no_std

Callers rarely see a Workspace, because a thread-local pool supplies one transparently. The common gemm call allocates at most once per thread, and reuses that buffer for every later call.

The pool is also re-entrancy-safe. Nested rayon can re-enter a GEMM on a thread already inside one. For example, a worker might work-steal another GEMM while blocked in its own for_each, or a batch-parallel worker might run an element inline. In that case the pool’s RefCell is already borrowed. So with_thread_pool hands out a fresh scratch workspace for that one call, instead of panicking. Packing buffers hold no result state between calls, so this fallback is invisible. Only that one call skips the buffer reuse.

For explicit control there is the *_with tier. Every entry has a variant (gemm_with, gemm_packed_b_with, and so on) that threads a caller-owned Workspace through. From the second sufficiently large call on, this gives zero heap allocation. It is the tool for hot loops of small products and for latency-sensitive code, and Workspace::with_capacity avoids even the first-call allocation spike.

Without std there is no thread-local storage, so with_thread_pool simply builds a fresh workspace per call. Because parallel requires std, there are no threads to re-enter either. A caller who wants reuse on such a build holds its own Workspace and uses *_with. no_std and WebAssembly recommends this same pattern.

Parallel Execution

Parallelism in gemmkit lives in one small layer, gemmkit/src/parallel.rs (layer L2), with 2 jobs. It decides how many workers a problem deserves, and hands which work to each of them.

Both decisions are deliberately conservative, because more threads are not free. The layer’s design starts from one observation: the wrong worker count loses more performance than the wrong schedule does. Both decisions are also shaped so that the numerical result never depends on either one.

The user-facing surface is a single enum. It is Parallelism::Serial, or Parallelism::Rayon(n). Rayon(0), the default, means auto.

Workload-aware worker resolution

Parallelism::resolve turns the request into an actual partition count. The request is only one input. The workload is the other.

First comes the total-work serial gate. Below GEMMKIT_PARALLEL_THRESHOLD (default 48*48*256 for m*n*k), everything stays serial, before the resolver even samples the core count. Forking rayon for a product that takes only microseconds would cost more than the product itself. The gate precedes the request check, so even an explicit Rayon(n) stays serial below it.

Above the gate, the resolver honors an explicit count, capped by the core count and by the number of available jobs. So Rayon(huge) can neither over-subscribe the machine nor over-allocate per-worker pack regions. Only the auto path is heuristic. This keeps forced widths exact for tests and for scaling diagnostics.

The auto path is work-based. It divides the total work m*n*k by GEMMKIT_PAR_MNK_PER_WORKER (default 2_000_000, the per-worker floor below which fork/join overhead outweighs the gain). This division gives the worker count, floored at 1 and capped by the core count and the job count.

The path is work-based rather than dimension-based, because the optimal worker count tracks total flops, not linear size. A small cube runs fastest serial. A mid-size cube wants only a few workers. A large cube wants every hardware thread the machine has. No single stride on a linear dimension can fit that spread.

Scaling to full width at mid sizes depends on one more thing: avoiding redundant per-worker packing of the same A panel. Packing and Workspaces covers the LHS in-place gate that keeps that redundancy from happening. The GEMMKIT_PAR_MNK_PER_WORKER knob is the escape hatch for a machine whose per-worker floor differs from the compiled default.

Bandwidth-bound shapes get a different rule entirely. A gemv or gevv does O(1) arithmetic per byte, so the compute ramp’s logic does not transfer. resolve_bandwidth gates on bytes touched instead.

Below a cache-derived byte floor, the matrix fits one core’s private cache, and that core saturates it alone. Splitting then only adds fork/join and shared-cache contention, with no bandwidth to gain. gemv_parallel_floor_bytes (in cache.rs) derives this floor from the topology. On parts with an L3, the floor is the per-core private L2. On parts without one, it is a fraction of the full shared cluster L2. GEMMKIT_GEMV_PARALLEL_BYTES overrides the floor directly.

Above the floor, the matrix spills to the shared L3, whose bandwidth a single core cannot saturate on its own. So the auto count steps straight to a wider width for those bytes, instead of ramping up to it. That width climbs a ladder built from the exact-fit pool tiers described below. The smallest tier sits at the floor, and the count climbs one tier for each GEMMKIT_GEMV_TIER_STEP (default auto, 8) factor of touched bytes. The ladder stops at the largest tier rather than the full machine width. A gemv saturates its bandwidth well before the machine runs out of cores, and more workers past that point can cost more than they gain.

gemmkit builds the rungs from the pool tiers themselves, not from a separate set of fractions. So an auto gemv width always has an exact-fit pool waiting for it, and it never pays the slack tax those tiers exist to remove. GEMMKIT_GEMV_THREAD_CAP replaces the whole ladder with one flat width, for a deployment that wants to fix it.

There is no ramp between rungs, because a bandwidth-bound scaling curve dips at a handful of workers. Fork/join and contention costs are already paid there, and aggregate bandwidth has not yet arrived. Any ramp through that dip would lose to both of its endpoints. So the rule stays simple: serial below the floor, a tier width above it, and nothing in between.

Batched GEMM has its own resolver, resolve_batch, which chooses between 3 plans. Serial runs every element in turn on the calling thread. BatchParallel(n) gives each worker whole, cache-hot GEMMs to run, paying one fork/join for the whole batch. Because no element ever splits, this plan stays bit-identical across worker counts. SequentialInternal instead loops the batch on the calling thread, giving each large, DRAM-bound element the full engine parallelism in turn.

resolve_batch gates the SequentialInternal split to m, n > 1 shapes, whose routes are worker-count independent. A gemv-shaped element always stays whole on one worker instead. Special Paths covers the routing, and Batched GEMM covers the API.

Demand-driven work distribution

Given a worker count, the driver does not build a nested task tree. For each column block and depth slice, it flattens the inner work into a flat 1-D job list: n_mc row blocks times n_nt column tiles. Job q decodes to (ic_idx, jt) = (q / n_nt, q % n_nt). Workers pull contiguous chunks from a shared, lock-free cursor until it is empty:

#![allow(unused)]
fn main() {
// gemmkit/src/parallel.rs
impl JobCursor {
    /// Atomically claim the next `[start, end)` chunk, or `None` once the job space
    /// is exhausted
    #[inline]
    pub(crate) fn next_chunk(&self) -> Option<(usize, usize)> {
        let start = self.next.fetch_add(self.grain, Ordering::Relaxed);
        if start >= self.n_jobs {
            None
        } else {
            Some((start, (start + self.grain).min(self.n_jobs)))
        }
    }
}
}

Each claim costs one fetch_add. There are no locks and no per-job queues. Demand-driven pulling is what makes heterogeneous cores work well. On a big.LITTLE part, a P-core that finishes chunks faster simply pulls proportionally more of them. A static n_jobs / n_threads split would instead leave every core waiting on the slowest one. The same property also absorbs OS noise and frequency differences on homogeneous machines.

The chunk grain balances 2 costs. Too coarse a grain leaves the tail of the job list idling workers at the join. Too fine a grain lets the atomic claims start to show, along with, on the packed-LHS path, re-packing at chunk edges.

The general grain oversamples the worker count. job_grain targets GEMMKIT_PARALLEL_OVERSAMPLE chunks per worker (default 8), so each worker expects several pulls, and imbalance self-corrects.

The packed-LHS path is special. Its natural chunk is a whole row block (n_nt consecutive jobs). A worker packs the block’s A panel once, and reuses it across all the block’s column tiles. That yields only n_mc chunks. So when the row-block count is small, packed_block_grain splits each block into power-of-two column sub-chunks. It splits until there are about GEMMKIT_PACKED_OVERSAMPLE * n_threads chunks (default target 2), and only by divisors of n_nt. So a chunk never straddles a row-block boundary and re-packs A mid-chunk. Splitting harder than this target re-packs too often and makes performance worse.

2 parallel phases run before the compute region in each depth slice, and their boundaries are the only barriers in the driver.

When B packs, workers pull nr-wide column panels from their own cursor. The fork/join of that region is the write-before-read barrier the compute region depends on, because packed B is the one buffer shared non-disjointly across workers.

The shared-A pre-pass packs each row block’s panel once into a shared slot, under the same discipline. It opens above a size gate, and it also opens from 16 workers up regardless of size. At that width, every extra worker is another redundant copy of each panel it touches, so deduplicating pays off even on mid-size problems.

Everything else is disjoint by construction. Workers write only their own output tiles and their own pack regions. That invariant is what lets the Ptr shim declare the captured raw pointers Send + Sync, in one audited place.

Size-class thread pools

Rayon’s fork/join cost does not scale with the problem. It scales with the pool’s idle slack: the gap between the threads a pool owns and the workers a given call actually engages. Forking k workers into a w-wide pool wakes w threads, not k. The w - k threads that get no work still pay their share of the barrier and the OS-level wake/park round trip.

The full-width global pool is the worst case for a small parallel GEMM. A mid-size product often wants only a fraction of the machine’s threads. Forking it into the full-width global pool then wastes most of that width as pure slack tax on every call.

gemmkit’s answer is a small set of private, persistent pools, each sized to exactly one of the worker counts the auto path actually asks for. Below full machine width, it keeps up to GEMMKIT_POOL_CLASSES halving tiers (default 2 on x86_64, 1 on aarch64, clamped to 3, 0 elsewhere). The tiers are a half-width pool, a quarter-width pool, and so on, each halving the machine’s physical width again. Each tier pool builds lazily, on first use, at a small one-time cost, and it never rebuilds afterward. The tiers are a fixed halving of the machine width, not a value tuned per shape.

The auto path snaps its worker count exactly onto a tier, so it carries zero slack by construction. It stays on its largest tier until the total work m*n*k clears GEMMKIT_FULL_WIDTH_MNK (default auto, arch-split: 110_000_000 on x86_64, 14_000_000 on aarch64). Past that point, the extra full-width workers finally pay for the fork/join they add, and full machine width takes over.

3 rules keep this mechanism from ever fighting the caller’s own scheduling.

1. A call already running inside a rayon pool, whether the caller’s own install or a nested gemmkit call, never redirects to a tier pool. The ambient pool always wins, exactly as it did before tier pools existed.
2. An explicit Rayon(n) keeps its exact semantics of precisely n workers. It only picks the smallest tier pool that fits n, instead of forking into the global pool. The worker count stays unaffected. Only the pool it forks into changes.
3. Threaded wasm keeps its own dedicated pool, described below, untouched by any of this. Tier pools are a native, non-wasm concern.

The wasm story

On wasm32-wasip1 there are no threads to spawn, and rayon would trap if it tried. The compile-time constant RAYON_USABLE captures whether the target can run workers at all. On a wasm build without the threading opt-in, every resolver returns 1, and for_each_worker runs the plain serial loop. parallel degrades gracefully instead of trapping.

The opt-in is the wasm_threads feature, which targets wasm32-wasip1-threads or a browser with SharedArrayBuffer. Because available_parallelism is unsupported on wasm, rayon’s global pool would otherwise silently size itself to 1 thread. So gemmkit builds its own pool instead, sized by the GEMMKIT_WASM_THREADS knob (default 8), and installs it around the worker loop. The deployer states the width, and everything else stays unchanged. no_std and WebAssembly has the details of the wasm builds.

The reproducibility contract, assembled

gemmkit’s reproducibility contract is simple. For a fixed machine and a fixed configuration, the engine produces reproducible output. This is not a promise of bitwise-identical results across different configurations, and the worker count counts as part of that configuration.

Today, gemmkit also holds a stronger property, as an engineering fact rather than a separate promise. Changing only the worker count, with everything else held fixed, still produces bitwise-identical output. 3 mechanisms make that true. Here is how they fit together.

First, the numerics do not depend on the worker count. Blocking derives kc and nc from the cache model alone, never from the thread count. The depth slices always run in the same fixed pc order. So every output element’s floating-point reduction takes the same shape, whether one worker drains the cursor or many do. The one dimension a worker count can move is mc. A wide count shrinks it, through the driver’s parallel job-depth floor, so the flat job list stays several chunks deep per worker. The job list itself is therefore not always identical, chunk for chunk, across widths. But mc always stays an mr multiple. So the microtile set (every mr-aligned row offset plus the one m-tail tile) is the same under any split. With kc and the pc order untouched, no result bit moves. The worker count changes how work is grouped and partitioned. It never changes what each tile computes.

Second, no reduction ever splits across workers. Within a depth slice, one worker computes an output tile’s whole update, so a chunk is always a set of whole tiles. The depth slices themselves run sequentially, since the pc loop is never parallel. beta applies on the first slice, and later slices accumulate. So every output element’s floating-point reduction runs in one fixed order, determined by the blocking alone. The special paths keep the same discipline. gemv partitions output rows on register-panel boundaries, so each row’s SIMD/scalar split does not depend on the partition. That makes gemv bit-identical across worker counts outright. Batched plans either keep whole elements on one worker, or split only shapes whose routes are already worker-count independent.

Third, packed bytes do not depend on who packs them. pack_panels is a pure reorder, and both of its branches write identical bytes. A panel packed by 1 worker, by the shared-A pre-pass, or by prepack_rhs earlier, is therefore the same sequence of bytes. See Packing and Workspaces. The kernel cannot observe who staged its input.

Which worker computes a given tile genuinely varies from run to run, since the cursor hands chunks to whoever asks first. By the 3 mechanisms above, though, nothing numerical depends on that choice. This bitwise agreement across worker counts is a property of gemmkit’s design today, not a wider guarantee. It leaves room for a tolerance-held kernel, such as the bf16 dot path, to reshape its accumulation without breaking the actual contract. Parallelism in Practice covers how to choose worker counts in practice, and what the contract means for testing.

Special Paths

The register-tiling driver described in Life of a GEMM Call is built around 1 assumption. It assumes there is enough work per output tile to amortize packing, blocking, and a full MR x NR register accumulator. Some shapes break that assumption badly. A matrix-vector product has no tile reuse at all. A k = 4 product finishes before the pack pays for itself. An 8 x 8 x 100000 contraction would spend most of the driver’s effort multiplying zero padding.

Layer L6 (gemmkit/src/special/) reroutes exactly these shapes to dedicated kernels. Every reroute holds 3 properties. It hides behind the same public entries, so gemm and its siblings never expose which path ran. It is threshold-gated through tuning knobs, so a mis-calibrated gate can move or turn off without a recompile. It also keeps the library’s reproducibility contract: the same machine and the same configuration give a reproducible result. Most of the routes below go further and stay bit-identical across worker counts for a fixed shape. gemv is the exception. Its own section below explains why.

The gates live at the top of each per-type dispatch entry, in a fixed order. This is run_typed from gemmkit/src/dispatch/float.rs, trimmed for clarity:

#![allow(unused)]
fn main() {
// gemmkit/src/dispatch/float.rs (run_typed, trimmed)
if (t.n == 1 || t.m == 1) && core::cmp::min(t.m, t.n) <= tuning::gemv_threshold() {
    gemv::run_typed_epi::<T, S, Identity>(/* user frame, before orientation */);
    return;
}
orient_transpose(&mut t);
if small_mn_eligible(&t) || small_mn_pack_eligible(&t) {
    small_mn::run_epi::<T, S, Identity>(/* horizontal dot kernel */);
    return;
}
if t.k <= tuning::small_k_threshold() {
    small_k::run::<FloatGemm<T>, S, MR_REG, NR>(/* one depth panel, in place */);
    return;
}
driver::run::<FloatGemm<T>, S, MR_REG, NR>(/* the general blocked driver */);
}

gemv fires before orientation normalization, in the user’s coordinate frame. The other gates run on the oriented problem. Each special path also exists in a fused-epilogue form, so a gemm_fused call takes the same route its plain twin would. That contract is the subject of Epilogue Fusion.

gemv: the memory-bound edge

A shape with m == 1 or n == 1 (gemmkit/src/special/gemv.rs) does 2k flops per output element, against k matrix elements read once. This makes it memory-bound. The whole design question is how to cut DRAM traffic, not how to schedule FMAs. Both orientations reduce to 1 core routine. It views the matrix, transposed when m == 1, as a rows x k block times a k-vector.

The gate checks shape, not size. min(m, n) is 1 for any gemv shape, so the GEMMKIT_GEMV_THRESHOLD knob it compares against works as an on/off switch rather than a size cap. Set it to 0 to force gemv shapes onto the general driver instead, which still produces a correct result.

Parallelism follows a bandwidth model, not the compute ramp. Parallelism::resolve_bandwidth stays serial below a cache-derived byte floor. Below that floor, the matrix fits inside 1 core’s private cache, and that core saturates it alone. Above the floor, the resolver steps straight to a width sized for the bytes touched. That width climbs a ladder over the exact-fit pool tiers. It tops out at half the logical cores, because a handful of workers is the worst point on a bandwidth scaling curve.

row_sweep partitions output rows across workers in panels whose grain is a multiple of the SIMD width. Each row’s full k-reduction stays inside 1 worker, and no worker ever merges a partial result from another. This keeps gemv reproducible at a fixed worker count, the same bar the rest of the engine holds to. Splitting the rows across more or fewer workers does not, by itself, promise bitwise agreement across those different worker counts. small-k, small-mn, and batched, below, go further than that.

Inside a worker’s row range, the code picks 1 of 4 strategies based on layout. A column-major matrix takes the axpy shape. It offers 2 variants that are deliberately bit-identical and differ only in memory traffic.

The register-blocked output form holds a panel of output rows in SIMD registers across the whole k-sweep. It reads the matrix and the output exactly once each. The plain column-outer form re-reads the output every few columns, but it streams the matrix as 1 contiguous read. The choice between them depends on output cache residency, computed in output_register_block. The route picks register-blocking under 2 conditions. First, the output (rows * sizeof) must outgrow a fraction of the last-level cache, so the plain form’s re-reads would otherwise reach DRAM. Second, k must sit at or below GEMMKIT_K_STREAM_MAX (default 32). Past that k, the register-blocked form’s many concurrent column streams start to thrash the hardware prefetcher. Both variants run the same ascending-k fused accumulation per element, and the same per-row SIMD-versus-scalar split. Switching between them never changes a single output bit. It only changes speed.

A row-major matrix takes a dot form instead. It register-blocks rows in groups of 4 to overlap FMA latency chains, while each row still runs the shared, fixed-order dot_contiguous reduction. Fully strided operands fall back to a scalar loop.

One shape fits 2 of those classifications at once, and the way that tie breaks decides whether anything vectorizes at all. The axpy forms vectorize over output rows, holding lanes of them in a register. The dot form vectorizes over k instead. A matrix of a single row has both its strides equal to 1, so column-major and row-major describe the same bytes equally well. This is the pure dot product, m == n == 1.

Handing that shape to the axpy form leaves its vector loop unreachable, since while i + lanes <= e never holds when e == 1. The whole reduction then falls onto the scalar remainder. axpy_yields_to_dot avoids this. It gives the sweep to the dot form whenever the row count is short of 1 SIMD register, and the dot form’s own strides hold. The dot form also carries a wider accumulator tree, so it is the more accurate of the two. This choice matters in practice. Column-major adapter libraries such as nalgebra and faer describe a row vector exactly this way. A caller reaching for a dot product would otherwise land on the slow classification by default.

Whether the rows split across workers is a separate decision from which strategy computes them. For a column-major matrix the answer is usually no. The output-row axis is that matrix’s inner, fastest-varying memory axis. Cutting it hands every worker a strided walk over the entire matrix, where each worker consumes only its own slice of every column. The serial route instead makes 1 sequential pass: row_sweep short-circuits to a single body(0, rows) call, with no blocking at all. That pass already runs near the achievable single-stream rate. Below some row count, extra workers have little to win and a great deal of sequentiality to lose.

GEMMKIT_GEMV_AXPY_PAR_MIN_ROWS holds that row-count floor. Below it, the axpy split stays serial regardless of the requested worker count. This floor deliberately exempts 2 routes. A row-major matrix gives each worker whole k-contiguous rows, so its stream stays sequential even when split, and splitting it keeps winning at every size. The mixed-precision twin’s widening axpy is compute-bound enough to scale on the same column-major stream, so it keeps splitting worthwhile too.

The mixed-precision twin, run_mixed (feature half), serves f16/bf16 gemv. It uses the same row partition and the same reproducibility argument as the float routine above. Every load widens to f32 through the KernelSimd<N, N, f32, N> seam, and the reduction runs in f32. The result rounds to the narrow type exactly once, at the store.

A single asymmetry follows from that single-rounding rule. The mixed axpy always register-blocks the output. The plain column-outer form would re-read and re-write the narrow output every column group. That would round it once per group, instead of once per element.

The mixed fused gemv does not route here at all, on purpose. The float fused gemv fuses by re-reading each stored output and mapping it in place. That is bit-exact only because the float output is the accumulator. A narrow output has already rounded once at the store, so reading it back and mapping it again would round it twice. Instead of threading the epilogue through every widening store, the mixed fused entry keeps gemv shapes on the general driver. The driver already applies the epilogue in f32, before its single narrowing step (gemmkit/src/dispatch/mixed.rs).

small-k: one depth panel, nothing to amortize

At a tiny k (gemmkit/src/special/small_k.rs), the whole product is a single depth panel. The driver’s cache-blocking model, its workspace carving, and above all its A/B packing would all be pure setup. Every packed element would be read only once anyway. This route instead computes C <- alpha*A*B + beta*C directly over the family’s microkernel, with kc = k. It reads A and B in place: no packing, no blocking, no workspace traffic. It still inherits the family’s widen and rounding semantics for free, because it stays generic over KernelFamily.

The gate is k <= GEMMKIT_SMALL_K_THRESHOLD. Its default splits by architecture: 16 on x86, 8 on aarch64. The narrower NEON microkernel tile packs cheaply enough that the driver wins sooner there, hence the lower default.

The in-place read needs 3 preconditions. When any of them fails, the route defers to the driver instead, and the result stays correct, just differently scheduled. First, the microkernel needs unit-stride LHS rows, so A must be column-major (rsa == 1). Second, a FORCE_PACK_* family, such as complex, transforms the data into a planar form while packing, so it cannot be read in place by construction. Third, k past a hard SMALL_K_MAX = 32 would overflow the one stack buffer this route uses: a zero-padded panel for the bottom, partial row-tile. That panel still needs packing, because the microkernel always loads a full mr rows.

The route partitions work over output tiles, one full k-pass per tile, with 1 worker per tile. The bandwidth model caps the worker count itself, since the m*n output write dominates at a small k. Every tile is a complete reduction owned by a single worker. So the result stays bit-identical across worker counts, a stronger property than gemv holds to.

small-mn: horizontal dots for tiny outputs

Both m and n can sit far below the microtile while k stays long (gemmkit/src/special/small_mn.rs). There, the driver would pad the tiny row and column tiles up to full MR x NR microtiles, and compute mostly padding. This route instead computes each output element as 1 horizontal SIMD dot, C[i,j] = alpha*<A[i,:], B[:,j]> + beta*C[i,j], streaming along the contraction. The output is register-blocked into 4 x 4 tiles of accumulators, so 16 independent FMA chains stay in flight across the k-sweep. Each A-row and B-column loads only once per tile. This is the same latency-hiding trick as gemv’s dot form, generalized to a small grid.

The dims gate is m, n <= GEMMKIT_SMALL_MN_DIM (default 16, 32 on aarch64), together with k above the small-k threshold. The cap is arch-split because the point where the driver’s padding starts to cost more than this route’s horizontal dots differs by machine. The 2 small-shape routes split the k axis between them this way.

The kernel needs both operands unit-stride along k. It wants A’s rows contiguous (csa == 1, row-major A) and B’s columns contiguous (rsb == 1, column-major B). This is the zero-copy tier. The 2 most common layouts each fail exactly 1 side: all-row-major fails B, and all-column-major fails A. For these, a second, mutually exclusive gate (k > GEMMKIT_SMALL_MN_PACK_MIN_K, default 16) engages a pre-pack tier. prepack_operands copies only the failing operand into a k-contiguous workspace scratch buffer, then runs the identical kernel over it with unit strides.

The copy touches m*k (or n*k) elements, against the m*n*k work of the dot itself. This is a small tax next to the horizontal kernel’s win, so a strided small-m,n shape still beats falling back to the driver’s padded microtiles. The scratch buffer rounds each line up to an odd number of cache lines (packed_line_stride). A natural stride of exactly k would map every packed line to the same L1 set whenever k is a power of 2. That collapses the benefit of the re-reads, so the odd-line rounding avoids it.

That tax is small in flops, but flops are the wrong unit for it. The copy does no arithmetic per byte it moves. The dots do about 2. The copy therefore has less to hide memory latency behind. On a long k it takes a far larger share of the time than of the work. On the calling thread it dominated the route.

The copy now runs across workers itself. It forks after its traffic clears the cache-derived byte floor that the bandwidth-bound routes share (GEMMKIT_GEMV_PARALLEL_BYTES). Below that floor the copy stays serial. It resolves its own worker count, apart from the tile sweep that follows. The 2 axes offer different amounts of parallelism. The MT x NT output grid caps the sweep, and a small m, n makes that grid tiny. The depth caps the copy, and a long k makes the depth large.

The copy splits the depth, never the lead axis. A contiguous t range lets each worker read whole depth lines, lead consecutive elements per step. A split of the few lead lines would instead take a single element per step from each line, across the whole operand.

Measured on the Zen5 reference machine (f32, auto width, column-major A): 8x8x524288 3.1x, 16x16x262144 2.0x, 4x4x1048576 1.8x, 8x8x2097152 1.7x, 16x16x1048576 1.1x. Footprint is the one trend that holds across these shapes. At a fixed m,n, the smaller packed operand gains more. The order across m,n at a fixed footprint is not monotonic. This page therefore claims no mechanism beyond the copy’s share of the serial time.

The pre-pack step is a pure reorder: the same values, in the same per-line order. So the packed route stays bit-identical to the eligible-layout route. Each worker writes every cell once, with the value a serial copy writes. A split of the copy across workers therefore moves no byte either. This route has 2 siblings. One is mixed (f16/bf16, widen to f32, round once per cell). The other is integer (i8 -> i32, wrapping, and therefore bit-exact against the driver). Both share the same tiling, the same pre-pack helper, and the same reproducibility argument.

batched: orchestration, not a kernel

Batched GEMM (gemmkit/src/special/batched.rs) is deliberately not a new kernel. Every element re-dispatches through the full single-GEMM engine, so a batch composes with the driver, gemv, small-k, and small-mn automatically. What this layer adds is a schedule, chosen once per call by Parallelism::resolve_batch:

• BatchParallel: chosen when there are enough elements to fill the workers. Each worker runs whole GEMMs serially and cache-hot, and the batch pays 1 fork/join instead of 1 per element. This is the model for the motivating workload: many small matrices.
• SequentialInternal: chosen for few, large, DRAM-bound elements. It loops the batch on 1 thread and gives each element the full engine parallelism in turn. On x86 the split engages once an element spills the per-core L2. aarch64 shares an L2 cache across a cluster and has high unified bandwidth. There, it engages once the per-batch-worker share, elem_bytes / batch, exceeds GEMMKIT_SEQ_INTERNAL_BYTES_PER_WORKER (default 128 KiB). This plan splits 1 element’s own work across workers, so it is gated to m, n > 1 shapes. The driver, small-k, and small-mn routes all reduce every output within 1 worker. So they agree bit-for-bit between serial and parallel, under the current, thread-independent blocking. gemv is held only to the base reproducibility promise, not to that bitwise agreement, so this plan excludes it.
• Serial: chosen below the total-work gate, or when no threads are usable.

Every element is independent, and the serial and batch-parallel plans never split one across workers. So a batch stays bit-identical across worker counts under those 2 plans, whatever the per-element route. The strided-batched entries, gemm_batched and gemm_batched_fused, thread 1 shared epilogue through the same skeleton and share a single schedule implementation. The plain and fused forms cannot drift apart.

For a batch whose elements differ in shape, the pointer-array form gemm_batched_ptr_unchecked takes a slice of GemmProblem descriptors instead. Each descriptor carries its own dimensions, strides, and pointers. This form uses the simpler resolve_batch_flat policy instead. It hands whole GEMMs to workers, and never splits a single element, since there is no uniform per-element residency to test. gemm_batched_slice is its safe, validated twin.

Batched GEMM and Small Shapes and GEMV cover the user-facing view of all this. Parallel Execution describes the worker-count machinery these paths depend on.

Epilogue Fusion

A GEMM output rarely leaves the routine raw. Inference layers add a bias and an activation. Quantized pipelines requantize the i32 accumulator down to a byte. Done naively, each of those is a second full pass over C: every element gets written to memory, evicted, read back, transformed, and written again.

The epilogue feature (gemmkit/src/kernel/epilogue.rs) fuses the transform into the microkernel’s store instead. It transforms the element in the register, or scratch slot, it already occupies, at the moment the microkernel would have stored it anyway. The second pass disappears. The saving is even larger for requantization. An unfused flow would have to materialize the entire m x n matrix in i32 before narrowing it.

The seam

The seam is the Epilogue trait. It threads through KernelFamily::microkernel_epi, so every family’s store site can apply it, and the driver never needs to know it exists:

#![allow(unused)]
fn main() {
// gemmkit/src/kernel/epilogue.rs (trimmed)
pub trait Epilogue<Fam: KernelFamily>: Copy + Send + Sync {
    /// true => every hook const-folds away; the kernel is bit-identical to non-fused
    const IS_IDENTITY: bool = false;
    /// true => apply_reg is implemented, enabling the fast vector store path
    const VECTOR: bool = false;
    /// true => apply_store is implemented (the Out != Acc requantize pattern)
    const VECTOR_STORE: bool = false;

    /// Scalar transform at absolute (row, col) in the oriented problem frame
    unsafe fn apply(&self, v: Fam::Acc, row: usize, col: usize) -> Fam::Out;
    /// Vector transform of LANES consecutive rows; MUST agree with apply bit-for-bit
    unsafe fn apply_reg<S>(&self, simd: S, v: ..., row: usize, col: usize) -> ...;
    /// Vector transform of a whole MR_REG x NR register tile; defaults to a loop over
    /// apply_reg. Overridden to hoist a runtime discriminant out of the unrolled pass
    unsafe fn apply_tile<S, const MR_REG: usize, const NR: usize>(&self, simd: S, acc: ..., row0: usize, col0: usize) -> ...;
    /// Vector store-transform from Acc scratch to Out; same bit-agreement contract
    unsafe fn apply_store<S>(&self, simd: S, src: *const Fam::Acc, dst: *mut Fam::Out, ...);
}
}

The whole design rests on 2 invariants. The first is zero-cost identity. Plain gemm passes the Identity epilogue, whose IS_IDENTITY = true makes every hook const-fold away. The monomorphized, non-fused kernel stays bit-identical to what it was before the seam existed. Fusion costs nothing when a caller does not use it.

The second is fire-once semantics. The driver hands the microkernel a last_k flag, and the epilogue applies only on the final depth panel. Earlier panels store raw Acc partials, exactly as the non-fused kernel would. A family with OUT_IS_ACC = false, such as the narrow f16/bf16 outputs, runs the whole contraction as a single kc = k panel, by construction. So last_k is structurally true there. The deep-K twin would break that single-panel guarantee, so it deliberately never engages on the fused path. The special paths fire once for free too, since each of their output elements is a single complete reduction with a single store.

The built-ins

The library ships 3 epilogues, each behind its own public entry point (feature epilogue). The requantize entries also need int8. See Fused Epilogues for the user-facing view:

FusedEpi is the runtime-composed bias-plus-activation epilogue. It adds a per-row or per-column bias (Bias::PerRow / Bias::PerCol), then applies Relu or LeakyRelu(slope). A single monomorphization covers every combination, so the fused kernel count does not multiply by the number of epilogue kinds. That only holds because the enum branches decode once per tile, in the apply_tile override, rather than once per accumulator.

That distinction is not a micro-optimization. The tile const generics unroll the kernel’s store pass. A per-register hook would replicate both match statements in every one of the tile’s accumulator slots. On a wide tile, the resulting branch web costs the compiler the accumulator tile itself. It stops keeping acc in registers. Instead, it writes every value through to the stack from inside the kc loop, where the epilogue does not even run. Hoisting the decode out of the loop, so it runs once per tile instead, avoids that spill entirely. tests/perf/fused.rs pins the ratio between the fused and the plain rate, so a future change cannot let this regression back in unnoticed.

FusedEpi backs gemm_fused and its whole constellation. This includes gemm_batched_fused, which shares 1 bias and 1 activation across the whole batch, the prepacked twins gemm_packed_b_fused and gemm_packed_a_fused, and the complex entry gemm_cplx_fused. The complex entry is bias-only, because an ordering-based activation has no mathematical definition on complex numbers. FusedEpi sets VECTOR = true. On the fast path, the bias add and the activation run as register operations, such as max(v, 0). Its NaN contract on the SIMD max/min is chosen so the vector and scalar forms agree exactly: both compute ReLU(NaN) = 0.

MapEpi is the escape hatch. gemm_map applies an arbitrary user closure, f(value, row, col) -> value, to each output element at its final value. It uses (row, col) in the user frame. The closure is a borrowed &dyn Fn + Sync, so there is 1 monomorphization per (type, ISA), not 1 per closure. It runs scalar, once per element, amortized by the O(k) flops behind each element. MapEpi supports f32/f64 only. A narrow type would have to round to N, apply the N-domain closure, then round again, which would break the bitwise contract described below.

KRequantize implements the quantized-inference store: C[r,c] = clamp(zp + round_ne(scale*(acc + bias)), LO, HI). It maps the i32 accumulator down to i8 (gemm_i8_requant, band [-128, 127]) or u8 (gemm_i8_requant_u8, band [0, 255], the ONNX QLinearMatMul convention). The scale can be per-tensor or per-row (RequantScale). The zero point joins in as an integer after the rounding step. An optional i32 bias joins in as an integer before the single f64 rounding step. The rounding itself is round-half-to-even, through a no_std-safe 2^52 trick (round_ne_f64). KRequantize has no alpha, since that folds into the scale, and no beta, since accumulating into an already-quantized C has no clear meaning.

The correctness contract

This contract is stated precisely because it is exactly what the epilogue tests pin down, bit by bit. It is composed of 3 ingredients:

1. Identical routing. A fused call routes every shape through the same kernel plain gemm would use. The general driver, gemv, small-k, and small-mn each exist in a fused form (see Special Paths). The fused dispatch entries mirror the plain gates one for one. No shape pays the driver’s overhead just because it asked for a bias. The one deliberate exception is the mixed f16/bf16 fused gemv, which stays on the driver for a rounding reason explained in Special Paths. Narrow types sit outside the bitwise contract below anyway.
2. An epilogue-independent engine. Blocking, scheduling, packing, and the accumulation order do not depend on which epilogue is threaded through. The epilogue only touches the store.
3. Bit-agreeing apply paths. A full column-major tile stores through the vector path, apply_tile, itself defaulting to apply_reg, or through apply_store. An edge or strided tile drains through scratch and the scalar apply instead. A single output matrix can freely mix the two. So the trait contract requires both paths to agree bit-for-bit under the same token. An apply_tile override inherits that obligation. It must leave exactly what apply_reg would have, element for element.

Together these 3 give the headline guarantee. For f32/f64, gemm_fused, gemm_map, and the batched and prepacked fused entries all equal gemm() followed by the same scalar map, bitwise, for every shape.

MapEpi shows how deliberate that guarantee is. It sets VECTOR = true. This is not to vectorize the closure, which it cannot do. It is so the kernel takes the same path selection plain gemm does. The fast path’s fused beta*C + alpha*AB store differs from the scalar path’s unfused arithmetic by 1 ULP for a general beta. A scratch-only epilogue would therefore hand the closure a value plain gemm never actually wrote. Instead, apply_reg drains the register to a stack buffer, and calls the same scalar apply once per lane. So f always sees exactly the bits plain gemm produced.

The documented exception is f16/bf16. The narrow blanket implementation applies the bias and the activation in f32, on the accumulator, before the single round-to-nearest-even narrowing to the output. That is deliberately more precise than gemm()-then-map, which would round to the narrow type, widen it back, and round again. So for narrow types the fused entries are not bitwise equal to gemm-then-map. The documentation states this plainly, rather than weaken the fused semantics to match the less precise alternative. Within 1 fused run, the vector and scalar paths still agree bit-for-bit, since both compute act(bias(v)) in f32 and round exactly once. Reproducibility across worker counts stays unchanged too.

KRequantize earns its vector path differently. The x86 tokens implement KernelSimd::requant_store: a vectorized widen-to-f64, scale, hardware round-to-nearest-even, clamp, and low-byte store. Its documentation carries a case-by-case proof that every lane equals the scalar clamp(zp + round_ne(scale*v), lo, hi). The i32 -> f64 and f32 -> f64 widenings are exact. The 2^52 trick agrees with the hardware rounding below 2^52, and saturation agrees above it. A NaN cannot occur, because the API validates that every scale is finite and positive. A per-row scale varies per lane, so that case takes the per-lane scalar map instead. Non-x86 tokens keep REQUANT_VECTOR = false and use the scalar map throughout. An in-module conformance sweep, the requant_store tests in gemmkit/src/simd.rs, checks this bit-equality on every capable token. The word “proven” in this contract is therefore enforced by a test, not just an aspiration.

One last corner remains. When the A*B term vanishes, because k == 0 or alpha == 0, the fused entries still owe C <- act(beta*C + bias). That degenerate map runs element-wise in the user frame. fused_degenerate, in gemmkit/src/dispatch/float.rs, handles this, with a narrow sibling that combines in f32 and narrows once. So even the no-op-product case honors the same semantics as the full kernel.

Extension Points

gemmkit’s variation points are traits, const generics, and typed function pointers stored in OnceLock slots. There are no macros and no transmute. That discipline exists for 1 reason. The library expects 4 kinds of growth: a new instruction set, a new element type, a new dot-product instruction, and a new fused transform. Each one should land as additive code, with a short, checkable list of touch points. Each one should also leave the driver, the packing routines, and the blocking model untouched.

This page turns those 4 recipes into walkthroughs. It is written for someone extending the crate itself. The seams stay public enough that the most important one, driving the generic driver with your own kernel family, also works from outside the crate. A test proves it.

A new ISA backend

An ISA backend is a zero-sized token, plus a set of vocabulary implementations. The wasm simd128 backend (gemmkit/src/simd/wasm.rs) is the most recent complete example. It is worth reading end to end, since it is only 1 file plus a handful of dispatch lines.

The token’s only inherent behavior is Simd::vectorize, the #[target_feature] trampoline. Runtime CPU detection cannot pair with a fixed #[target_feature] attribute on the generic kernel. So every kernel invocation runs inside a tiny, annotated function. The #[inline(always)] primitives fold into that function, so every intrinsic lands in feature-enabled codegen:

#![allow(unused)]
fn main() {
// gemmkit/src/simd/wasm.rs
impl Simd for Simd128 {
    #[inline(always)]
    unsafe fn vectorize<R>(self, f: impl FnOnce() -> R) -> R {
        #[target_feature(enable = "simd128")]
        fn inner<R>(f: impl FnOnce() -> R) -> R {
            f()
        }
        inner(f)
    }
}
}

The checklist:

1. The token. Add a Copy + Send + Sync + 'static zero-sized struct in a new gemmkit/src/simd/ module. cfg-gate it to its architecture, and give it the vectorize trampoline shown above.

2. SimdOps<T> implementations. Add 1 for each element type the ISA accelerates. Each one needs a register type and LANES. It also needs the primitive vocabulary: load, store, splat, mul, add, mul_add, fnma, reduce_sum, plus max/min if the fused float epilogue should vectorize. The vocabulary stays deliberately thick, so the microkernel can stay 1 generic function. You implement primitives, never a kernel.

   Honor the documented contracts here. The simd128 implementation uses f32x4_pmax, not f32x4_max, because the trait’s max requires a NaN in a to yield b. This is the ReLU(NaN) = 0 agreement between the vector and scalar epilogues. It also passes the operands reversed, as f32x4_pmax(b, a), because pmax(x, y) computes x < y ? y : x. The natural argument order would return NaN for a NaN a, and -0.0 for max(-0.0, +0.0), the opposite of the contract in both cases. It also spells mul_add as an unfused mul followed by add, because wasm has no hardware FMA. The relaxed-SIMD alternative is nondeterministic by specification, which would break reproducibility.

3. Tile geometry. Pick (MR_REG, NR) for each type, and encode it as the const generics of the per-ISA wrapper functions in the dispatch modules. This is the only per-(type, ISA) knob. Budget the registers explicitly. simd128 runs 2x4 for f32: 8 accumulators, 2 LHS registers, and 1 RHS register, for 11 live v128 values. LLVM’s wasm backend starts to spill past about 16 live vectors, which is why simd128 stays at that width. NEON runs 4x4 instead, leaving registers spare on purpose.

4. A Dispatched descriptor, and 1 arm per select_* ladder. The memoized selection ladders live in gemmkit/src/dispatch/. They are select_f32/select_f64 for float, select_f16/select_bf16 for mixed, select_i8 for int, and select_c32/select_c64 for complex, plus the map-epilogue selectors. Each ladder arm bundles the plain, prepacked, and fused entry points together with the tile geometry. So adding the ISA costs 1 descriptor constant and 1 match arm per type it accelerates.

5. A GEMMKIT_REQUIRE_ISA name. Add a ForcedIsa variant and its parse string in gemmkit/src/dispatch/isa.rs. The current values are scalar, fma, avx512f, avx512vnni, avx512bf16, neon, simd128, and auto. Follow the fail-loudly rule. If the pinned ISA is unsupported, dispatch must panic rather than fall back. This way, a CI job that means to exercise your kernel cannot silently pass on a different one instead.

6. Tests ride along mostly for free. tests/simd_conformance.rs constructs tokens directly and checks every primitive against scalar references. An env_isa_* pin binary, plus a CI job, makes the dispatch route itself testable too (see Testing and Verification).

You should not need to touch driver.rs, any kernel family, pack.rs, or cache.rs. The simd128 backend changed none of them.

A new element type

Element types vary along 2 small traits (see Scalars and Kernel Families). Scalar (gemmkit/src/scalar.rs) declares only the identity constants and the accumulator type Acc. Choosing Acc is the single most consequential decision here, since it fixes the rounding story. f16 chose Acc = f32. i8 chose Acc = i32, which makes integer GEMM exact. KernelFamily (gemmkit/src/kernel.rs) bundles everything else that distinguishes the operation: the Lhs/Rhs/Acc/Out types, the pack layout, and the microkernel.

Often a new family is not needed at all. If the type is a narrow input over an existing accumulator, implement the KernelSimd<L, R, A, O> widen/narrow seam on the capable tokens instead. That means widening loads, plus 1 narrowing store. Then reuse the generic microkernel, the way MixedGemm<f16> and MixedGemm<bf16> do. The homogeneous case is a blanket implementation, and the mixed implementations cannot overlap it. A genuinely new operation shape, such as the planar complex kernel or the requantizing integer families, gets its own KernelFamily instead.

Wiring the new type into the public API means adding a dispatch module under gemmkit/src/dispatch/, with its own OnceLock slot per type. Feature detection runs once, the winning monomorphized entry points get cached, and every later call becomes 1 indirect call. Copy dispatch/mixed.rs as a pattern for 2 types, gemv/small-mn/small-k reroutes, and a dot-kernel selection wrinkle. Copy dispatch/int.rs instead for a heterogeneous task type.

The open/closed property here is not folklore. gemmkit/tests/open_closed.rs enforces it directly. That test defines NaiveFloat, an independently written family with its own packing and a plain scalar microkernel, built using only public items. It then drives the unchanged public driver::run with NaiveFloat, and checks the result against an f64 reference. If a driver change ever breaks the family seam, that test fails to compile. It also works as the template to start a new family from.

A dot-product instruction

Instructions such as vpdpbusd and vdpbf16ps fold several depth steps into 1 operation. That reshapes the accumulation rounding, so they must never arrive as a clever override of the portable tile loop. The seam is split to keep that distinction clear:

• The family declares DEPTH_MULTIPLE = Q (greater than 1), and packs through pack_kgroup_panels (gemmkit/src/pack.rs). That function interleaves Q consecutive depth steps contiguously per lane. The driver rounds panel depths up to Q, and keeps k-groups from straddling slice boundaries.
• The capable token overrides KernelSimd::dot_accumulate, and consumes whole instruction groups from those panels. The packed layout is a private contract between the family’s packers and the overriding token. Any signedness correction, such as VNNI’s +128 trick with its column-sum compensation, lives inside the override. So the accumulator holds the true sum when it returns.
• SimdOps::accumulate_tile overrides are reserved for scheduling changes that keep the rounding shape unchanged. 2 examples are an in-order core that needs explicit software pipelining, and a scalable-vector ISA whose length is not a compile-time constant. Its documentation is explicit that rounding-reshaping instructions are out of scope for this seam. Those should arrive as a new family with the dot seam instead. An accumulate_tile override must stay deterministic, and must round consistently with the edge path. The default implementation already saturates the FMA pipes on any wide out-of-order core, so prove an override earns its keep before adding one.

IntGemmVnni and Bf16DotGemm are the 2 worked examples. IntGemmVnni is bit-exact against the widen path, because integer arithmetic is associative. Bf16DotGemm is held to a tolerance instead, within the reproducibility contract. Dot Kernels and the Deep-K Twin covers both in depth.

A new fused transform

A fused transform is just 1 Epilogue implementation (gemmkit/src/kernel/epilogue.rs). The driver’s last_k plumbing, the zero-cost Identity default, and the routing through every special path all come for free (see Epilogue Fusion). The design work is choosing an application path, and honoring 1 hard rule. The vector and scalar paths must agree bit-for-bit. Full tiles take the vector path, edge and strided tiles take the scalar path, and a single output matrix can mix the two.

• A transform on Acc-typed values with a natural register form sets VECTOR = true and implements apply_reg. This is the FusedEpi pattern. Mind the NaN and signed-zero semantics here. LeakyRelu, for instance, is written as the identical max + slope*min composition in both forms.
• A transform that narrows Acc to a different Out sets VECTOR_STORE = true and implements apply_store. This is the KRequantize pattern. Argue the bit-equality case by case, and pin it with a conformance sweep.
• A transform with no profitable vector form keeps both flags false, and routes everything through scratch and the scalar apply. This stays correct for any tile shape. If the scalar value could differ by 1 ULP from the fast path’s fused store, though, borrow MapEpi’s trick instead. Set VECTOR = true, and implement apply_reg as a drain-to-stack-and-apply-per-lane. This way the transform always sees exactly the bits plain gemm would have stored.

Whatever path a new transform takes, add its own gemm-then-map equivalence test next to the existing ones in gemmkit/tests/epilogue/. That suite is where the bitwise contract actually gets enforced.

Testing and Verification

A library whose headline promises are “bit-identical here, tolerance there, reproducible everywhere” lives or dies by how precisely its tests pin those words down. gemmkit’s suites live in gemmkit/tests/. The first structural decision is what is not a test. The performance harnesses are measurement tools. They never gate CI.

tests/perf/ is the exhaustive internal investigation suite. It runs #[ignore] benchmarks over a median-of-9 harness, serialized behind a shared lock, because each one saturates every core. Someone runs it by hand when a change needs numbers. gemmkit/benches/gemm_bench.rs is the curated public cargo bench surface instead. It holds criterion benchmarks in 5 headline groups, sgemm, dtypes, gemv, prepacked, and batched, meant for --save-baseline regression tracking against the gemm crate and matrixmultiply. Neither suite can fail a merge. A performance assertion on a shared CI runner would mostly assert noise.

Correctness, properties, conformance, fuzzing

The correctness suite (tests/correctness/) sweeps shapes, layouts, and alpha/beta combinations against an independent f64 reference GEMM. The reference, and the accuracy machinery around it, live once in tests/oracle_common/. This includes the element traits, the deterministic fills, the f64 reference itself, and the relative-Frobenius accuracy gates for each element type. Both the correctness and the property harnesses include this module with #[path], so there is exactly 1 oracle to trust.

On top of the oracle sweeps sit several more checks. Cross-checks against the external gemm crate catch shared-blind-spot bugs the in-repo reference cannot, since it is an independent implementation. Parallel bit-identity tests cover the routes where the library promises that agreement. Per-ISA kernel runs go through the generic driver. The safe API’s exact panic wording, easy to underrate, gets held by #[should_panic(expected = ...)] substrings, such as "A.cols" and "aliases itself". So a validation message cannot quietly degrade into a less useful one.

The property tests generalize these sweeps, and all 3 drive proptest over shapes, strides, and knob values. tests/props_api.rs covers oracle accuracy, run-to-run bit determinism, serial-equals-parallel agreement, the beta == 0 overwrite semantics, broadcast strides, batching, and panic guarantees. tests/props_packed.rs covers prepacked-versus-plain bit-identity in the general regime, and tolerance on the documented tiny/gemv exception set. tests/props_knobs.rs covers behavior under randomized knob settings.

One layer down, tests/simd_conformance.rs checks the L0 vocabulary itself. It constructs every ISA token the host supports directly, bypassing dispatch. It then compares each SimdOps primitive, the homogeneous KernelSimd blanket, and the portable fma_bvec fallback, lane-by-lane, against scalar references. This is where primitives the product kernels rarely touch, such as integer reduce_sum, fnma, and the widen seam, get exercised at all. In-module sweeps, such as the requant_store bit-equality tests in gemmkit/src/simd.rs, do the same for the vector requantize contract. The suite has no proptest dependency, so it also runs on wasm and conformance-tests the compile-time simd128 token.

Fuzzing lives in gemmkit/fuzz/, a nightly-only cargo-fuzz sub-workspace with its own workspace root, excluded from the stable build. It holds 6 libFuzzer targets. fuzz_gemm builds valid-by-construction problems and checks them differentially against naive references, so any panic there is a library bug. fuzz_batched does the same for batched calls. fuzz_prepack and fuzz_prepack_i8 round-trip through the prepack APIs, with the i8 one gated bit-exactly. fuzz_api_validation throws adversarial geometry at the checked entries. There, a documented "gemmkit:" panic counts as an accepted outcome, and anything else counts as a validation gap. fuzz_knobs sets every process-global tuning knob to an adversarial value before each run. This is the target that mechanically finds arithmetic-overflow classes in the blocking model.

Isolation discipline

A naive test layout turns racy around 2 kinds of global state, and the suite’s structure exists to work around both.

Tuning knobs are process-global atomics. Every test that mutates one lives in its own dedicated binary. tests/tuning.rs holds the setters, tests/env.rs holds environment-variable resolution, and there is also tests/props_knobs.rs and tests/deep_k_narrow.rs, which toggles GEMMKIT_DEEP_KC_BYTES to force each deep-k route. tests/env.rs holds exactly 1 test, so its environment access is single-threaded by construction. A separate binary is a separate process, and it cannot race another binary’s knob state.

Within one binary, though, libtest still runs tests concurrently. So every knob-touching test there serializes under a per-binary KNOB_LOCK mutex, and restores whatever it changed before it releases that lock. The property-test binary adds an RAII guard on top, one that survives proptest’s internal catch_unwind.

GEMMKIT_REQUIRE_ISA is stickier still, because dispatch memoizes it once per process. So there is 1 pin binary per value: tests/env_isa_avx512f.rs, _vnni, _bf16, _scalar, _neon, and _wasm, plus env_isa_garbage.rs, which asserts the unknown-value panic. Each binary routes every test through a shared Once in tests/env_isa_common/. That Once performs the single set_var call before any dispatch resolves. Because every test in one binary pins the same value, it does not matter which test wins the race to run that Once. The write deliberately overrides an inherited GEMMKIT_REQUIRE_ISA. That override is what lets the SDE-pinned CI jobs below run these same binaries and still exercise the real, per-ISA routes.

Miri rounds out the memory-safety story where fuzzing’s sanitizers stop. CI runs the scalar-path correctness suite (miri_scalar_path) and the complex negative-stride unchecked entry under Miri. Miri interprets the actual unsafe pointer arithmetic of the pack and microkernel paths directly. A cfg(miri) detour exists only where Miri cannot interpret a hardware conversion. It never exists to skip logic.

The CI matrix

.github/workflows/ci.yml turns this pinning machinery into coverage of kernels the runners do not physically have:

SDE emulation runs far slower than native execution. This is where GEMMKIT_FAST_TEST earns its keep. It is a test-suite-only switch, implemented once in tests/fast_test_common/ and included from the harnesses. The library itself never reads it. The switch shrinks the deterministic dimension and coefficient sweeps down to 1 representative per redundant combination, while still visiting every branch and path class. The SDE jobs set it, alongside PROPTEST_CASES=16. Native jobs keep the full sweeps instead. Keeping this switch out of the library proper means a test-convenience flag can never change shipping behavior.

The net effect ties this chapter together. Every claim the earlier pages made is held by a test you can point to, in a binary whose isolation rules make the result trustworthy. This includes the bit-identity guarantees in Special Paths, the gemm-then-map equivalence in Epilogue Fusion, and the open/closed property in Extension Points. It also includes each pinned kernel’s correctness on hardware the project does not own.