CPU 后端

主线位置：算子主线在 CPU 上的终点。Storage::matmul 匹配到 (Cpu, Cpu)（见 Storage 与后端 trait），调用 CpuStorage::matmul，它只有一行：MatMul(bmnk).map(self, lhs_l, rhs, rhs_l)。本页讲这一行背后的两层：Map2::map 按 dtype 分派，MatMul::f 把数交给 gemm crate 或 Apple Accelerate 的 sgemm；顺带讲逐元素运算和线程池。

flowchart TB
  s["Storage matmul"] --> c["CpuStorage matmul"]
  c --> m2["Map2 map 按 dtype"]
  m2 --> f["MatMul f T"]
  f --> q{"feature"}
  q -->|默认| g["gemm crate"]
  q -->|accelerate| a["sgemm dgemm"]
  q -->|mkl| k["MKL gemm"]
  g --> pool["rayon 当前线程池"]
  g --> v["Vec T"]
  a --> v
  v --> w["CpuStorage F32 等"]

代码在 candle-core/src/cpu_backend/（mod.rs 3327 行，utils.rs 363 行，conv2d.rs）和 candle-core/src/cpu/（SIMD 小函数，neon.rs、avx.rs、simd128.rs），Accelerate 的 FFI 声明在 accelerate.rs，线程池在 utils.rs。

本页的实验 crate（M2 Max，8 个性能核 + 4 个能效核；macOS 27、rustc 1.99.0；全部 release）：exp-candle/backend-ops（candle-core 0.11.0，features metal accelerate）和 exp-candle/backend-noaccel（不开任何 feature）。

CpuStorage：每种 dtype 一个 Vec

CpuStorage 是一个枚举，变体就是 dtype：

#[derive(Debug, Clone)]
pub enum CpuStorage {
    U8(Vec<u8>),
    U32(Vec<u32>),
    I16(Vec<i16>),
    I32(Vec<i32>),
    I64(Vec<i64>),
    BF16(Vec<bf16>),
    F16(Vec<f16>),
    F32(Vec<f32>),
    F64(Vec<f64>),
    F8E4M3(Vec<F8E4M3>),
    // Dummy types that store raw bytes
    F6E2M3(Vec<u8>),
    F6E3M2(Vec<u8>),
    F4(Vec<u8>),
    F8E8M0(Vec<u8>),
}

数据就是一个普通的 Vec<T>，元素类型由变体决定，不需要像 Metal 那样另存一个 dtype 字段。最后四个变体存原始字节，几乎所有运算对它们都返回 UnsupportedDTypeForOp（为什么有这四个 dtype 见 DType 与 Device）。Clone 是派生的，try_clone 和 to_cpu_storage 都是 Ok(self.clone())，即整块复制。

每个运算都返回一个新的 CpuStorage，也就是一个新分配的 Vec。没有缓冲池，内存复用交给系统分配器。

Map1 / Map2：dtype 分派写一次

一个 14 个变体的枚举，两个输入就是 196 种组合。cpu_backend/utils.rs 用几个辅助 trait 把这个 match 只写一遍：

pub trait Map2 {
    const OP: &'static str;
    fn f<T: WithDType>(&self, v1: &[T], l1: &Layout, v2: &[T], l2: &Layout) -> Result<Vec<T>>;

    fn map(&self, v1: &C, l1: &Layout, v2: &C, l2: &Layout) -> Result<C> {
        match (v1, v2) {
            (C::U8(v1), C::U8(v2)) => Ok(C::U8(self.f(v1, l1, v2, l2)?)),
            (C::U32(v1), C::U32(v2)) => Ok(C::U32(self.f(v1, l1, v2, l2)?)),
            // ...
            (C::F32(v1), C::F32(v2)) => Ok(C::F32(self.f(v1, l1, v2, l2)?)),
            (C::F64(v1), C::F64(v2)) => Ok(C::F64(self.f(v1, l1, v2, l2)?)),
            (C::F8E4M3(v1), C::F8E4M3(v2)) => Ok(C::F8E4M3(self.f(v1, l1, v2, l2)?)),
            _ => Err(Error::DTypeMismatchBinaryOp {
                lhs: v1.dtype(),
                rhs: v2.dtype(),
                op: Self::OP,
            }
            .bt()),
        }
    }
}

实现者只写泛型的 f<T>，拿到的是已经解开的 &[T] 和各自的 Layout；map 负责解包、调用、再包回同一个变体。一个算子是一个结构体，参数放在字段里：MatMul((b, m, n, k))、Affine(mul, add)、IndexSelect { ids, .. }。同一个文件里还有几种变体：

trait	签名差别	用户
`Map1`	一个输入，输出同 dtype	`Affine`、`ReduceSum`、`IndexSelect`、`Gather`、池化、上采样、`Im2Col`
`Map1Any`	输出 dtype 由实现决定	`ReduceIndex`（argmin / argmax 输出 `u32`）
`Map2`	两个输入，同 dtype	`MatMul`、`WCond`、`IndexAdd`、各种卷积
`Map2InPlace`	第一个输入可写	`Scatter`
`Map2U8`	输出 `u8`	`Cmp`（比较）

dtype 不同在这里会报 DTypeMismatchBinaryOp，但多数情况到不了这一步，Storage 那一层的 same_dtype 已经拦下了。

逐元素的 unary_impl / binary_impl 没有用 Map：它们要的不是一个泛型 f<T>，而是 B::f32、B::f32_vec 这些按 dtype 分开命名的函数（见 Storage 与后端 trait 里的 UnaryOpT），所以 CpuStorage::binary_impl 自己手写了一个按变体的 match，每支调 binary_map_vec(lhs_l, rhs_l, lhs, rhs, B::f32, B::f32_vec, B::f32_scalar_vec)。

按布局遍历：逐元素运算

binary_map_vec 先用 NdIter 把两个布局里能合并的维度合并，得到一个「最内层切片」的长度和两边在这一层的步长，再按步长分四种情况：

    for [lhs_off, rhs_off] in nd_iter {
        match (inner_ls, inner_rs) {
            (1, 1) => f_vec(
                &lhs[lhs_off..lhs_off + inner_size],
                &rhs[rhs_off..rhs_off + inner_size],
                &mut ys_to_set[dst_off..dst_off + inner_size],
            ),
            (1, 0) => {
                let r = rhs[rhs_off];
                f_scalar_vec(
                    r,
                    &lhs[lhs_off..lhs_off + inner_size],
                    &mut ys_to_set[dst_off..dst_off + inner_size],
                );
            }
            (0, 1) => {
                let l = lhs[lhs_off];
                for i in 0..inner_size {
                    ys_to_set[dst_off + i] = f(l, rhs[rhs_off + i]);
                }
            }
            _ => {
                for i in 0..inner_size {
                    ys_to_set[dst_off + i] =
                        f(lhs[lhs_off + i * inner_ls], rhs[rhs_off + i * inner_rs]);
                }
            }
        }
        dst_off += inner_size;
    }

两边都连续时调切片函数 f_vec，右边是广播的标量时调 f_scalar_vec，其余逐元素调标量函数 f。输出 Vec 用 with_capacity 分配后直接写未初始化的部分，最后 set_len，省掉一遍清零。一元运算的 unary_map_vec 换成按 Layout::strided_blocks() 分块，道理一样；StridedBlocks 的三种形态见布局与视图。

这两个函数里没有 rayon。CPU 上的逐元素运算、to_dtype、affine、普通的归约都是单线程的；用到 rayon 的是 matmul（经 gemm）、卷积（into_par_iter）、candle-nn 里的 softmax（par_chunks），以及下面要讲的量化矩阵乘和 CPU flash attention。

cpu/ 里的 SIMD

f_vec 的默认实现是一个 iter().zip().for_each()，靠编译器自动向量化。手写的 SIMD 在 candle-core/src/cpu/，入口是 VecOps trait：vec_add、scalar_add、vec_dot、vec_reduce_sum、vec_reduce_max、vec_reduce_min，每个都有标量的默认实现，按类型覆盖。f32 的 vec_add 是三选一：

    fn vec_add(lhs: &[Self], rhs: &[Self], res: &mut [Self]) {
        #[cfg(feature = "mkl")]
        crate::mkl::vs_add(lhs, rhs, res);
        #[cfg(all(feature = "accelerate", not(feature = "mkl")))]
        crate::accelerate::vs_add(lhs, rhs, res);
        #[cfg(not(any(feature = "mkl", feature = "accelerate")))]
        lhs.iter()
            .zip(rhs)
            .zip(res)
            .for_each(|((&a, &b), y)| *y = a + b)
    }

架构相关的部分按 target_feature 选：aarch64 上是 neon.rs 的 CurrentCpu，Unit = float32x4_t、STEP = 32、EPR = 4，一次循环处理 32 个 f32、用 8 个累加寄存器（ARR = STEP / EPR），最后两两相加归约。vec_dot_f32 和 vec_sum 建在它上面，f16 / bf16 的 vec_add、vec_dot 也在这里（先加宽成 f32 再算）。

这些函数的用户不是 matmul，而是行内的点积和归约：ReduceSum 在「沿最后几维求和、输入连续」时逐行调 vec_reduce_sum；candle-nn 的 softmax 用 vec_reduce_max 和 vec_reduce_sum（见自定义算子）；CPU flash attention 用 vec_dot。

matmul：gemm 或 Accelerate

Map2 for MatMul 有三个 f，用互斥的 cfg 选一个。默认的那个（第 1373 行）调 gemm crate（0.19.0）：

        match T::DTYPE {
            DType::F16 | DType::F32 | DType::F64 => {}
            _ => Err(Error::UnsupportedDTypeForOp(T::DTYPE, "matmul").bt())?,
        }
        // ...
        let (a_skip, b_skip) = self.ab_skip(lhs_l, rhs_l)?;
        let c_skip: usize = m * n;
        // ...
        let mut dst = vec![T::zero(); b * m * n];
        let num_threads = crate::utils::get_num_threads();
        let parallelism = if num_threads > 1 {
            Parallelism::Rayon(num_threads)
        } else {
            Parallelism::None
        };
        let (b, m, n, k) = if b_skip == 0 && a_skip == m * k {
            // a_skip and c_skip should be updated but step is always 0 so
            // it wouldn't matter.
            (1, b * m, n, k)
        } else if a_skip == 0 && b_skip == n * k {
            (1, m, b * n, k)
        } else {
            (b, m, n, k)
        };
        for step in 0..b {
            // ...
            unsafe {
                gemm(
                    /* m: usize = */ m,
                    // ...
                    parallelism,
                )
            }
        }

几个要点：

dtype 是白名单：F16、F32、F64。BF16 不在里面，CPU 上 bf16 张量做 matmul 直接报错（实测见下）。这个白名单是 #813（「Proper error on unsupported dtypes when using gemm.」）加的，在那之前只有 BF16 被单独拒绝，整数 dtype 会一路传进 gemm。
gemm 吃的是行步长和列步长，所以转置过的输入不用先 contiguous()；ab_skip 只要求批维能用一个步长描述，否则报 MatMulUnexpectedStriding。
一边的批步长为 0（广播）时，把批维并进 m 或 n，一次 gemm 算完；否则按批循环，每批一次 gemm。并进 m（右边广播）没问题；并进 n（左边广播、a_skip == 0）时输出按 (m, b·n) 排，和 [b, m, n] 的布局对不上，0.11.0 上这种输入会静默算错（broadcast_as 出来的左边再 matmul 就会走到这里），开了 accelerate 的 BLAS 路径和 Metal 不受影响。main 上的 #3960 修了它，还没有发布。
线程数来自 get_num_threads()，交给 gemm 的 Parallelism::Rayon(n)；gemm 在当前 rayon 线程池里用 into_par_iter 开 n 个任务。

开了 accelerate 时用另一个 f，调 Accelerate 框架的 BLAS sgemm（Fortran 接口 sgemm_）：

                        crate::accelerate::sgemm(
                            transa, transb, /* m= */ n as i32, /* n= */ m as i32,
                            /* k= */ k as i32, /* alpha= */ 1., /* a= */ a,
                            /* lda= */ lda, /* b= */ b, /* ldb= */ ldb,
                            /* beta= */ 0., /* c= */ c, /* ldc= */ n as i32,
                        )

BLAS 按列主序存矩阵，candle 按行主序。行主序的 C = A·B 在列主序眼里就是 Cᵀ = Bᵀ·Aᵀ，所以这里把 a 和 b 对调（a 是 rhs 的指针、b 是 lhs 的），m 和 n 也对调，算出来的列主序 Cᵀ 正好是行主序的 C。转置标志 transa / transb 由最后两维的步长判断，只认「连续」和「刚好转置」两种，别的步长报 non-contiguous。这条路只支持 F32、F64；F16 写的是 crate::bail!("the accelerate backend does not support f16 matmul")，BF16 和整数报 UnsupportedDTypeForOp。mkl 的那个 f 结构相同，多支持 F16（hgemm）。

Accelerate 也接管了一部分逐元素运算：上一页讲过 unary_op! 在开 accelerate 时把 Exp、Log、Sin、Cos、Tanh、Sqr、Sqrt 的 f32_vec 换成 vvexpf 这类 vForce 函数（vs_exp），Add 的 vec_add 换成 vDSP 的 vDSP_vadd。

实测：开与不开 accelerate

cpu_bench 在两个 crate 里是同一份源码，每项先跑一次预热，再计 9 次取中位数；矩阵是 randn 生成的 f32 方阵，逐元素运算的输入是 4194304 个 f32：

for n in [256usize, 1024, 2048] {
  let a = Tensor::randn(0f32, 1., (n, n), &dev)?;
  let b = Tensor::randn(0f32, 1., (n, n), &dev)?;
  let ms = median_ms(9, || a.matmul(&b))?;
  // ...
}
println!("exp  f32 [{n}]: median {:.3} ms", median_ms(9, || x.exp())?);

accelerate=true get_num_threads=8
matmul f32 256x256x256: median 0.032 ms (1042 GFLOP/s)
matmul f32 1024x1024x1024: median 0.867 ms (2477 GFLOP/s)
matmul f32 2048x2048x2048: median 7.843 ms (2190 GFLOP/s)
matmul f16 1024: Err(the accelerate backend does not support f16 matmul)
exp  f32 [4194304]: median 1.765 ms
add  f32 [4194304]: median 0.733 ms
mul  f32 [4194304]: median 0.624 ms
accelerate=false get_num_threads=8
matmul f32 256x256x256: median 0.174 ms (193 GFLOP/s)
matmul f32 1024x1024x1024: median 3.256 ms (659 GFLOP/s)
matmul f32 2048x2048x2048: median 28.961 ms (593 GFLOP/s)
matmul f16 1024: median 1.603 ms
exp  f32 [4194304]: median 6.408 ms
add  f32 [4194304]: median 0.580 ms
mul  f32 [4194304]: median 0.570 ms

项	accelerate	gemm（不开）	倍数
f32 matmul 1024	0.87 ms	3.3 到 4.5 ms	约 3.7 到 5
f32 matmul 2048	7.8 ms	29 ms	3.7
f16 matmul 1024	报错	1.6 到 2.5 ms	—
`exp`	1.8 ms	6.4 ms	3.6
`add`	0.73 到 0.88 ms	0.55 到 0.70 ms	Accelerate 略慢
`mul`（两边都没有覆盖）	0.53 到 0.62 ms	0.54 到 0.65 ms	相同

范围取自下一节五组线程设置里的各次运行。matmul 和 exp 是 Accelerate 真正快的地方；add 换成 vDSP_vadd 之后反而比自动向量化的 Rust 循环慢一点，4M 个元素的加法本来就受内存带宽限制，换函数省不下什么。gemm 路径上 f16 比 f32 快，是因为数据量减半；Accelerate 路径上 f16 干脆不能乘。

同一组 2×2 矩阵在 backend-ops 和 backend-noaccel 里按 dtype 各乘一次：

# backend-ops（accelerate）
cpu matmul F32: Ok([2, 2])
cpu matmul F64: Ok([2, 2])
cpu matmul F16: Err(the accelerate backend does not support f16 matmul)
cpu matmul BF16: Err(unsupported dtype BF16 for op matmul)
cpu matmul U32: Err(unsupported dtype U32 for op matmul)

# backend-noaccel（gemm）
F16: Ok([[2.0, 2.0], [2.0, 2.0]])
BF16: Err(unsupported dtype BF16 for op matmul)

也就是说，同一份用 f16 权重的代码，开 accelerate 之后在 CPU 上会从「能跑」变成报错；bf16 在 CPU 上两种构建都乘不了，Metal 上则三种浮点类型都可以（见 Metal 后端）。

线程：三个池

utils.rs 里和线程有关的有三样东西，读的环境变量也不一样：

pub fn get_num_threads() -> usize {
    // Respond to the same environment variable as rayon.
    rayon_num_threads()
}

名字	线程数	用户
`get_num_threads()`	`RAYON_NUM_THREADS`，没设时 `default_num_threads()`：macOS 上是性能核数（`hw.perflevel0.logicalcpu`，M2 Max 为 8），其他系统是物理核数	gemm matmul 的 `Parallelism::Rayon(n)`
`candle_pool()`	同上	只在 `Device::with_context` 里面生效
`barrier_pool()`	`CANDLE_NUM_THREADS`，没设时同上；主线程算一个，另起 n−1 个	量化矩阵乘（`k_quants.rs`）、CPU flash attention

此外还有 rayon 自己的全局池：卷积和 softmax 的 par_* 在没有 with_context 时跑在它上面，大小是 RAYON_NUM_THREADS 或逻辑核数（这台机器 12）。

Device::with_context 在 CPU 上就是 candle_pool().install(f)：

    /// Run `f` with device specific context.
    ///
    /// On CPU this installs candle's private rayon thread pool for the
    /// duration of `f`, keeping worker threads warm across the many short
    /// parallel sections in a model forward pass. Currently noop for other backends.
    pub fn with_context<F, R>(&self, f: F) -> R

这个池的工作线程在 macOS 上启动时调 pthread_set_qos_class_self_np(QOS_CLASS_USER_INTERACTIVE, 0)，让系统优先把它们放在性能核上。install 期间所有「当前池」的 rayon 调用（包括 gemm 的 into_par_iter）都进这个池。with_context 是可选的：candle 内部没有任何地方调用它，candle-examples 里只有 quantized-qwen3 一个例子包了一层。

barrier_pool 是 BarrierPool，一个手写的常驻线程池：每个工作线程一个独占缓存行的 Slot，主线程按 CPU 簇（cluster）只通知每簇的一个 leader，leader 再分给簇内其他线程，完成时反方向汇总；空闲时先自旋再 park。它为「一次前向里成百上千个很短的并行段」设计，rayon 的任务窃取在这种粒度上开销太大。

📝NOTE

冷知识：默认线程数在 0.11.0 里改了三次

0.10.x 里 get_num_threads() 没设环境变量时返回 num_cpus::get()，即逻辑核数，M2 Max 上是 12。PR #3584（「[CPU] Shared threadpool + less contention」，提交说明第一条是「Add Candle specific rayon threadpool. Use less cores by default (for now)」）加了 candle_pool 和 with_context，默认值改成性能核数，超过 4 个时再减半，注释写着「When the CPU backend is optimized further so as not to trash shared caches etc the amount of threads can be increased」，M2 Max 上变成 4。PR #3627 加了 BarrierPool，并让 CANDLE_NUM_THREADS 优先于 RAYON_NUM_THREADS。六天后的 PR #3643（「[CPU] Quantized perf improvements」，其中一条是「Have candle barrier pool respond to different env var than rayon」）删掉了减半，把两个环境变量拆开：get_num_threads 只看 RAYON_NUM_THREADS，barrier_pool 只看 CANDLE_NUM_THREADS。三个 PR 都首次进入 0.11.0，所以从 0.10 升到 0.11，gemm 的默认线程数从 12 变成 8。

实测：线程数

同一个 cpu_bench，用 RAYON_NUM_THREADS 设成 1、4、8、12（它同时决定 get_num_threads() 和 rayon 全局池的大小），1024 的 matmul 和 exp：

== cpu_bench_noaccel RAYON_NUM_THREADS=1
matmul f32 1024x1024x1024: median 21.769 ms (99 GFLOP/s)
exp  f32 [4194304]: median 6.499 ms
== cpu_bench_noaccel RAYON_NUM_THREADS=4
matmul f32 1024x1024x1024: median 8.061 ms (266 GFLOP/s)
exp  f32 [4194304]: median 6.480 ms
== cpu_bench_noaccel RAYON_NUM_THREADS=8
matmul f32 1024x1024x1024: median 4.479 ms (479 GFLOP/s)
exp  f32 [4194304]: median 6.513 ms
== cpu_bench_noaccel RAYON_NUM_THREADS=12
matmul f32 1024x1024x1024: median 3.636 ms (591 GFLOP/s)
exp  f32 [4194304]: median 6.618 ms

`RAYON_NUM_THREADS`	gemm 1024	Accelerate 1024	`exp`（gemm 构建）	`exp`（Accelerate 构建）
1	21.8 ms	0.90 ms	6.50 ms	1.80 ms
4	8.1 ms	0.91 ms	6.48 ms	1.84 ms
8（默认）	3.3 到 4.5 ms	0.86 ms	6.51 ms	1.89 ms
12	3.6 到 4.3 ms	0.88 ms	6.62 ms	1.81 ms

gemm 路径随线程数扩展，1 到 8 个线程快了约 5 倍；8 和 12 在三次重复里互有胜负（8：3.25、3.73、4.18 ms；12：3.78、3.88、4.35 ms），4 个能效核帮不上忙。
Accelerate 路径完全不受 RAYON_NUM_THREADS 影响：sgemm 用的是 Accelerate 自己的并行（和 Apple 的矩阵协处理器），candle 的线程设置管不到它。
exp 在两种构建里都不随线程数变化，它是单线程的。

with_ctx 程序比较「直接调」和「包在 Device::Cpu.with_context 里调」，各 201 次取中位数，不开 accelerate，两次运行：

available_parallelism=12 candle get_num_threads=8
n=64: plain median 0.0076 ms, with_context median 0.0075 ms
n=256: plain median 0.2198 ms, with_context median 0.1650 ms
n=1024: plain median 3.3249 ms, with_context median 3.1906 ms
available_parallelism=12 candle get_num_threads=8
n=64: plain median 0.0075 ms, with_context median 0.0075 ms
n=256: plain median 0.1828 ms, with_context median 0.0943 ms
n=1024: plain median 3.3073 ms, with_context median 4.1276 ms

64×64 太小，gemm 低于并行阈值，两边一样；256 的 matmul 在 with_context 里快了 25% 到 50%；1024 两次结论相反，落在噪声里。with_context 帮的是中等大小、调用频繁的并行段，这和它的文档注释一致；对单个大 matmul 看不出差别。

❗IMPORTANT

动脑筋：不开 accelerate，把 RAYON_NUM_THREADS 从 8 改成 1，4194304 个 f32 的 x.exp() 会慢几倍？

用上面的 cpu_bench_noaccel 跑两次，只看 exp 那一行。

动脑筋答案

== cpu_bench_noaccel RAYON_NUM_THREADS=1
exp  f32 [4194304]: median 6.499 ms
== cpu_bench_noaccel RAYON_NUM_THREADS=8
exp  f32 [4194304]: median 6.513 ms

一倍也不慢：6.50 ms 对 6.51 ms。Tensor::exp → CpuStorage::unary_impl::<Exp> → unary_map_vec，整条路径上没有 rayon，线程设置不起作用。
同一组运行里 1024 的 matmul 从 4.5 ms 变成 21.8 ms，因为 gemm 拿到的是 Parallelism::None。
想让 exp 变快，靠的不是线程而是开 accelerate（1.8 ms，单线程的 vvexpf）。

小结

CpuStorage 是按 dtype 分变体的 Vec<T> 枚举，每个运算分配一个新 Vec，没有缓冲池。
Map1 / Map2 等辅助 trait 把「解包 dtype、调泛型 f<T>、包回去」写一次，算子只写 f<T>；逐元素运算直接用 UnaryOpT / BinaryOpT 的按 dtype 函数。
逐元素运算单线程，按 NdIter / StridedBlocks 分出连续块，连续块调切片函数；cpu/ 里的 neon / avx 代码服务于点积和归约。
matmul 默认走 gemm（F16 / F32 / F64，按步长读，不必先连续），开 accelerate 走列主序的 sgemm（只有 F32 / F64）；bf16 在 CPU 上乘不了。实测 Accelerate 让 f32 matmul 快 3.7 倍、exp 快 3.6 倍，add 反而略慢。
gemm 的线程数是 RAYON_NUM_THREADS 或性能核数；量化矩阵乘和 flash attention 用另一个池、另一个环境变量 CANDLE_NUM_THREADS；Device::with_context 是可选的，对中等大小的 matmul 有帮助。