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.