CPU 后端The CPU backend

We covered ops saying WHAT matmul to do (L11), but how does it finally get computed, step by step, on a CPU with real machine instructions? This lesson drops to the lowest level and watches ggml's CPU backend (ggml/src/ggml-cpu/) take one dot product / matmul from "scalar, one multiply at a time" up to "SIMD, one instruction does a whole row", then split it across threads. This is the most hardcore stretch of the whole guide - but also the best place to see where performance actually comes from.

We use the quantized dot product as our thread - the hottest kernel in inference. It is the innermost loop of matmul, run countless times for every token the model generates; make it fast and the whole model is fast. This lesson reads real AVX2 SIMD code line by line; do not be scared - with the diagrams you will find its core idea is actually plain. You will gradually see that "hardcore" is hard not because of what any line computes, but because you must hold several threads in your head at once - how data is laid out, how instructions run in parallel, how the cache hits - and a diagram is exactly the tool that lets you see all those threads at a glance.

Roadmap: first the "scalar vs SIMD" difference (one at a time vs eight at once), then how the quantized dot product unpacks 4-bit weights and vectorizes them, and finally how multithreading splits one big matmul so many cores compute together.

From scalar to SIMD

First the plainest way. A dot product is "multiply matching positions and sum": sum += a[i]*b[i], looped n times. This is the kernel of the scalar reference ggml_vec_dot_f32 in ggml-cpu/vec.cpp - correct and readable, but slow: the CPU processes one number per cycle while its wide execution units sit mostly idle. This scalar version has an often-overlooked use too: it is the "correctness baseline" for every SIMD specialization. Any architecture's vector implementation must line up with it - bit for bit for the integer quantized kernels, within a tiny rounding error for float (multiple accumulators reorder the summation) - and anything off by more is a bug. So when reading low-level kernels, understand this slow-but-correct scalar version first, then read the fast SIMD version, and you hold a yardstick for "is it right" - a very practical reading order: grasp correct first, then grasp fast.

// scalar: one at a time (simplified from ggml-cpu/vec.cpp ggml_vec_dot_f32)
float sum = 0;
for (int i = 0; i < n; i++)
    sum += a[i] * b[i];                       // one multiply-add, repeated n times

// SIMD / AVX2: one instruction does 8 (vec.cpp GGML_SIMD path)
__m256 acc = _mm256_setzero_ps();             // 8-lane float accumulator
for (int i = 0; i < n; i += 8) {
    __m256 va = _mm256_loadu_ps(a + i);       // load 8 a's at once
    __m256 vb = _mm256_loadu_ps(b + i);       // load 8 b's at once
    acc = _mm256_fmadd_ps(va, vb, acc);       // acc += va*vb, 8 lanes at once
}
float sum = hsum_float_8(acc);               // horizontal sum of 8 lanes -> scalar

The SIMD (Single Instruction Multiple Data) on the right exists to squeeze out that idle compute. AVX2 offers 256-bit __m256 registers, one holding exactly 8 floats; one _mm256_fmadd_ps (fused multiply-add) instruction makes these 8 lanes each do one acc[i] += a[i]*b[i] simultaneously. The loop stride goes from 1 to 8, cutting instruction count by seven-eighths. ARM's NEON is 128-bit, 4 floats at a time - the same idea at half the width.

Seeing that "half the width" in real code is clearest - the same dot product, NEON just swaps 8 lanes for 4 and _mm256_* for v*_f32:

// NEON / ARM: one instruction does 4 (vec.cpp NEON path; macros in simd-mappings.h)
float32x4_t acc = vdupq_n_f32(0.0f);          // 4-lane float accumulator
for (int i = 0; i < n; i += 4) {
    float32x4_t va = vld1q_f32(a + i);        // load 4 a's at once
    float32x4_t vb = vld1q_f32(b + i);        // load 4 b's at once
    acc = vfmaq_f32(acc, va, vb);             // acc += va*vb, 4 lanes at once (FMA)
}
float sum = vaddvq_f32(acc);                  // horizontal sum of 4 lanes -> scalar

Map it line-by-line to AVX2: vdupq_n_f32(0) ~ _mm256_setzero_ps, vld1q_f32 ~ _mm256_loadu_ps, vfmaq_f32 ~ _mm256_fmadd_ps, vaddvq_f32 ~ hsum_float_8 - all the names change, the skeleton is identical. That is SIMD's portability: switching architecture only swaps one set of intrinsics; "load a row, multiply-add a row, finally horizontal-sum" stands unchanged.

Drawing this "8 lanes in parallel" is the clearest. Below we trace one SIMD dot product: 8 pairs multiply-add into 8 accumulators at once, and after a few rounds a single horizontal sum (hsum) folds the 8 accumulators into the final scalar.

How the quantized dot product works

In real inference, weights are quantized (L29), so the dot product must unpack before computing. Take the most common vec_dot_q4_0_q8_0: weights are Q4_0 - grouped 32-at-a-time into 16 bytes of 4-bit values (block_q4_0 = one fp16 d (scale) + qs[16]), activations are Q8_0 int8. So the inner of one quantized dot product is this chain: unpack 4-bit -> subtract 8 offset -> multiply int8 activation -> multiply back the scale -> accumulate. Below is its real AVX2 implementation, line by line.

// real AVX2 quantized dot product core (arch/x86/quants.c vec_dot_q4_0_q8_0)
__m256 acc = _mm256_setzero_ps();
for (; ib < nb; ++ib) {                               // loop over blocks (32 weights each)
    __m256  d  = _mm256_set1_ps(dx * dy);            // combine the two blocks' fp16 scale
    __m256i qx = bytes_from_nibbles_32(x[ib].qs);    // unpack: 16 bytes -> 32 values [0..15]
    qx = _mm256_sub_epi8(qx, _mm256_set1_epi8(8));   // offset to [-8..+7]
    __m256i qy = _mm256_loadu_si256((const __m256i*)y[ib].qs); // load 32 int8 activations
    __m256  q  = mul_sum_i8_pairs_float(qx, qy);     // int8 dot product -> float
    acc = _mm256_fmadd_ps(d, q, acc);                // FMA: acc += d * q
}
float sumf = hsum_float_8(acc);                      // horizontal sum of 8 lanes -> scalar

Line by line: acc is an 8-lane float accumulator; each loop iteration processes one block - bytes_from_nibbles_32 unpacks 16 bytes into 32 values in [0..15] (mask the low 4 bits with 0xF, shift for the high 4); subtract 8 to offset into [-8..+7] (4-bit quants are signed); _mm256_loadu_si256 loads 32 int8 activations at once; mul_sum_i8_pairs_float does the int8 dot product into 8 floats; finally _mm256_fmadd_ps multiplies by the combined scale and accumulates into acc. After all blocks, hsum_float_8 horizontally sums the 8 lanes into the final scalar. This "unpack + vectorized multiply-add" is exactly what lets a quantized model run on CPU. One easy-to-miss point: why are activations Q8_0 (int8) and weights Q4_0 (4-bit), different precisions on the two sides? Because their roles differ - weights are fixed and bulky, so squeezing to 4-bit saves the most memory; activations are live with a dynamic range, so int8 keeps them stable. And at the hardware level, int8 multiply-add has dedicated fast instructions (like _mm256_maddubs_epi16), even faster than a pure-float dot product. So "weights 4-bit, activations int8" both saves memory and runs fast, the mainstream of community quantization - and the very reason this kernel unpacks first, then does int8.

Pulling one block's flow out and freezing it makes each step clearer:

Multithreading: many cores together

SIMD squeezes a single core; multithreading puts all cores to work. One big matmul has many output rows, and each row's computation is independent (row 5 does not need row 3's result), so it is naturally parallel: split the rows evenly among N threads, each computes its batch, then merge. This "no dependencies, split however you like" structure is the ideal target for data parallelism. Hidden here is a plain but important test: whether you can parallelize comes down to "are there dependencies". Rows have no ordering between them, so you can split freely; once there is a dependency (a later step waits on an earlier result), parallelism needs synchronization and waiting, and the payoff drops immediately. When we look at the GPU in L32, you will find the very same test at work.

ggml's implementation is lightweight: each op ggml_compute_forward_* (ggml-cpu.c) gets params->ith (which thread am I) and params->nth (how many threads), and from them works out "which rows/blocks I own", computing on its own. The thread pool lives in ggml-threading, avoiding the cost of creating and destroying threads over and over. SIMD (8 at once within a core) x multithreading (many cores at once) stacked together is the entire source of the CPU backend's throughput - no magic, just spreading the same work out to be done as simultaneously as possible. Worth noting: not every op splits as nicely as matmul. Reduction ops like softmax and RMSNorm need to "see the whole row before computing", so splitting them needs care at the boundaries; elementwise ops (add, activations) split as freely as matmul. ggml writes per-op splitting logic exactly to handle these differences. Reading the source, you will see almost every ggml_compute_forward_* open by using ith/nth to work out its own range - that is what multithreading looks like landing at the op level.

Stacking the three: from dot product to matmul

We have seen SIMD (8 lanes within a core) and multithreading (many cores at once), and met tiling (reuse cached data, detailed in the fold below). But in one real matmul, the three are not independent - they nest layer upon layer and act together. See this relationship clearly and you truly understand "where the CPU backend's speed comes from".

Stacked, the effect is multiplicative: multithreading spreads work over all cores, tiling keeps each core from waiting on memory, and SIMD turns each core's "arm" into an 8-head one. Drop any layer and the whole thing slows - only multithreading without vectorizing leaves single cores crawling; only vectorizing without tiling leaves data shuttling endlessly between memory and cache. Because all three are on, ggml can run a multi-GB quantized model respectably on a pure-CPU machine with no GPU at all.

This also explains the tuning lore: thread count is usually best near the physical core count (more just fights for resources); and switching the build from no-SIMD to AVX2 or even AVX-512 often doubles speed at a stroke - now you know it is because the "inner" arm went from 1-at-a-time to 8 or 16. Low-level kernels look arcane, but the rule is concrete: whichever layer is not full, fill it; to find which layer, go back to L30's two rulers and measure.

Deep dive: cache and many architectures

Two final folds for two easy-to-overlook engineering details that truly make the CPU backend fast: what tiling does for the cache, and how one codebase fits all kinds of CPU.