逐元素、SIMD 与融合

主线位置：直接计算主线里 engine 上的逐元素运算（x.add(&y, &mut engine)），以及 traced 主线执行阶段的逐元素区域。前者从 CpuEngine 进 tenferro-internal-cpu-kernels（超越函数在 tenferro-cpu 的 analytic.rs），再进 strided-rs；后者由运行期把一串逐元素指令交给 tenferro-cpu-fused。两条路的输出都从 tenferro-cpu-basic 的缓冲池里取，用完再回去。engine 和令牌见 engine 与 threadteam 池，运行期怎么把图切成区域见编译和执行与缓存。

flowchart TB
  eng["CpuEngine add_read"] --> k["internal-cpu-kernels"]
  k -->|按 DType 匹配| typed["typed_add_with_pool T"]
  typed --> out["PooledUninitOutput"]
  typed --> er["erased_zip_into_uninit"]
  er -->|按 KernelDType 匹配| zm["zip 遍历"]
  zm -->|值得分派| team["team section"]
  an["exp tanh ..."] --> simd["elem-simd-math"]
  rt["traced 运行期区域"] --> fu["cpu-fused"]
  fu -->|不满足条件| fb["逐条执行"]
  fu --> sf["strided-fused"]
  sf --> spec["静态特化"]
  sf --> chunk["分块寄存器重放"]
  sf --> interp["逐元素解释器"]
  out --> pool["BufferPool"]

下面是 strided-rs 0.9（strided-kernel、strided-basic、strided-fused、strided-nn）和 elem-simd-math 0.1.2。实验在 exp-tenferro8/cpu-exp，release，M2 Max（8 个性能核 + 4 个能效核）、macOS 27、rustc 1.99.0，uptime 负载 2.2 到 2.8，时间取最小值。

一次 add 的路径

CpuEngine::add_read 铸一个令牌，调 add_read_with_pool，把缓冲池和令牌一起交下去。

第一次分派：tenferro 里按 DType。 两个输入都是拥有型 Tensor 时进 add_with_pool，对 (lhs.dtype(), rhs.dtype()) 做 match：六个同类型分支（F32、F64、I32、I64、C32、C64），加四个「实数标量配复数张量」分支，其余报错；有一个是视图时由 dispatch_read_presets! 展开出同样的组合。每个分支最后调同一个泛型函数，落到 zip_with_pool。它只允许 0 维操作数广播；一般的广播在上一层，TypedTensor::add 先用 engine 的 broadcast_in_dim_read 把形状不同的一侧展开。然后把输入转成 strided 的类型擦除描述符，交给 run_into：

fn run_into<O>(
    op: &'static str,
    dtype: DType,
    buffers: &mut BufferPool,
    shape: &[usize],
    write: impl FnOnce(&mut ErasedRawStridedUninitMut<'_>) -> Result<(), StridedError>,
) -> crate::Result<TypedTensor<O>>
where
    O: PoolScalar + KernelStorageElement,
{
    let strides = col_major_strides(shape)?;
    let mut out = PooledUninitOutput::<O>::new(buffers, shape)?;
    {
        let mut dest = ErasedRawStridedUninitMut::from_uninit_slice(
            out.as_uninit_slice_mut(),
            shape,
            &strides,
            0,
        )
        .map_err(|err| crate::Error::backend_source(op, err))?;
        write(&mut dest).map_err(|err| strided_error(op, dtype, err))?;
    }
    // SAFETY: a successful erased strided entry initializes every element of its dense column-major destination, and the destination descriptor is dropped above.
    unsafe { out.assume_init_recycled() }
}

输出永远是紧凑列主序，输入可以是任意步长的视图。注释说这是「the only place the elementwise family finalizes pooled output」，最后一行的 assume_init_recycled 下文再讲。

第二次分派：strided-rs 里按 KernelDType。 write 闭包调 strided-kernel 的 erased_zip_into_uninit，参数是运行期的 KernelDType 值、ErasedZipOp::Add 运算码和令牌。它核对 dtype、拒绝输入和输出重叠，再对 KernelDType 和运算码各 match 一次，进一个单态化的遍历。为什么在 tenferro 里拿到了 T 又擦掉、到 strided 里再匹配？擦除入口让「每种 dtype × 每种运算」的循环只在 strided-kernel 里单态化一次，tenferro 的几个 crate 只调编好的入口；运行时多出的是每次运算两次 match，不在元素循环里。

什么时候并行

0.7.1 的 strided 有一个固定门槛：超过 32768 个元素才切给线程池。strided-basic 0.9 换成一个代价模型，模块注释（cost.rs）的第一段是：「A kernel does not own a minimum length. It states how long its region takes on one thread (elements times the serial throughput of its class) and asks the placement token whether that pays」。每个内核属于四个 KernelClass 之一，每类在 aarch64-apple 上有一个每元素的串行耗时（注释说是在一颗 M5 Max 核上量的）：

令牌拿这个估计去比：串行时间不少于 3 倍「分派代价 × 分派次数」才并行（threadteam 的 PAYS；分派代价表见 engine 与 threadteam 池）。ew_sweep 程序在串行和 4、8 线程的 engine 上对不同长度的 f64 做 add（µs）：

add f64, min us:
  n     8192:  serial     1.62  team 4     2.08  team 8     2.12
  n    16384:  serial     3.29  team 4     3.88  team 8     4.00
  n    32768:  serial     6.08  team 4     6.58  team 8     6.71
  n    65536:  serial    11.50  team 4    12.00  team 8    12.17
  n   131072:  serial    22.25  team 4    13.88  team 8    15.67
  n   262144:  serial    44.42  team 4    21.83  team 8    14.46
  n  1048576:  serial   169.25  team 4   120.25  team 8    63.21

65536 个元素以下 team engine 只多出铸令牌的约 0.5 µs，说明没有分派；131072 起才并行。按 Map 的 0.1 ns 估计，65536 个元素约 6.6 µs、131072 个约 13 µs，门槛落在两者之间。实际的 add 比估计慢（65536 个要 11.5 µs，约 0.18 ns 每元素），所以门槛偏保守。

n    4096:  serial tanh   29.46  serial exp    7.71  serial mul    1.08  team 8 tanh   27.92  team 8 exp    7.50  team 8 mul    1.46  (min us)
n   16384:  serial tanh  109.29  serial exp   27.96  serial mul    2.96  team 8 tanh  109.83  team 8 exp   28.33  team 8 mul    3.58  (min us)
n   65536:  serial tanh  435.46  serial exp  109.67  serial mul   10.67  team 8 tanh  435.96  team 8 exp  110.46  team 8 mul   11.08  (min us)
n  262144:  serial tanh 1741.38  serial exp  436.67  serial mul   40.71  team 8 tanh  469.42  team 8 exp  120.38  team 8 mul   20.25  (min us)

超越函数走 SIMD 内核

exp、log、sin、cos、tanh、expm1、log1p 和 erf（f32、f64）在 simd_math::apply 里交给 elem-simd-math。这个库是 SLEEF 的 Rust 移植，基于 fearless_simd，运行时选指令集。模块注释给的理由是「Rust does not vectorize calls into libm, so those paths run at the speed of the scalar C library」。pow 故意不在其中（SLEEF 的 pow 比平台 libm 慢），sqrt / rsqrt 本来就会被 LLVM 向量化。

内核不在 tenferro 里决定怎么切块。map_unary 把它作为回调交给 map_into_uninit_with_kernel：strided 负责布局和并行，单线程时把整个缓冲区作为一段交给内核，并行时每个线程的每一块各交一次；内核拒绝的段（比如复数）退回逐元素闭包。这是 #2097 的做法，它的说明是让内核「now runs inside the parallel leaf instead of beside it」；之前 #2090 的第一版在 strided 旁边自己判断能不能走快路径，#2094 修了快路径把并行关掉的问题。三个都首次进入 v0.8.0。

上游的数字是 x86（AVX2）对 glibc 的 libm 量的。本机对 macOS 的 libm，ew_sweep 的后半段在 2^22 个元素上比较「普通的标量循环」和 tenferro（ns / 元素）：

analytic ops, n = 2^22, ns per element (min of 10):
  exp f64     std loop  2.211  serial  1.776  team 4  0.457  team 8  0.232
  tanh f64    std loop  3.507  serial  6.841  team 4  1.821  team 8  0.895
  sin f64     std loop  2.059  serial  2.436  team 4  0.641  team 8  0.323
  log f64     std loop  2.619  serial  3.057  team 4  0.807  team 8  0.404
  sqrt f64    std loop  0.307  serial  0.298  team 4  0.109  team 8  0.142
  exp f32     std loop  1.496  serial  0.569  team 4  0.149  team 8  0.274
  tanh f32    std loop  1.801  serial  2.903  team 4  0.751  team 8  0.386
  log f32     std loop  1.933  serial  1.364  team 4  0.362  team 8  0.189

• 单线程时，f64 只有 exp 比 libm 快（1.78 对 2.21）；tanh 慢了将近一倍（6.84 对 3.51），sin、log 也略慢。f32 的 exp、log 快 1.4 到 2.6 倍，tanh 慢 1.6 倍。
• 为了排除 tenferro 自己的开销，simd_direct 直接调 elem-simd-math 的切片函数：tanh f64: elem-simd-math 6.849 ns/elem, libm loop 3.490 ns/elem，exp 是 1.792 对 2.237，sin 是 2.430 对 2.098。差距来自内核本身。aarch64 上 f64 的向量只有两条通道，macOS 的 libm 又不慢，在这台机器上 SLEEF 式的 tanh 不划算。
• 多线程时这些都被并行掩盖了：8 线程的 tanh 是 0.9 ns，比单线程 libm 快 3.9 倍。

激活函数（sigmoid、silu、softplus、gelu、gelu_tanh）和 softmax 另有融合的 SIMD 内核（activation.rs、softmax.rs，后者调 strided-nn），一次遍历算完不物化中间结果，见 nn 运算。

输出：整块覆盖守卫与回池

PooledUninitOutput 是 run_into 用的守卫：

pub struct PooledUninitOutput<T: PoolScalar> {
    pool: BufferPool,
    shape: ShapeVec,
    strides: tenferro_tensor::StrideVec,
    data: Vec<std::mem::MaybeUninit<T>>,
    checkout: Option<UninitCheckoutToken>,
    byte_len: usize,
}

• new 从池里按「容量不小于需要的最小一档」取一块缓冲（take_best_fit），取不到就新分配，都不清零。写满之前调用方只能拿到 MaybeUninit<T> 的切片。
• 唯一的出口是 unsafe fn assume_init_recycled：调用方保证每个元素都写过，守卫用同一块分配构造 TypedTensor，并挂一个指向池的回收器（Weak<dyn HostBufferRecycler<T>>）。张量的最后一个所有者 drop 时，缓冲回到池里；池已经不在了（engine 被丢掉），就正常释放。
• 没写完就被丢掉的守卫（出错、或者故意放弃），Drop 把缓冲还给池：「leaves storage that is still the pool's, so it goes back for the next call rather than to the C allocator」。

BufferPool 现在是 Arc<SharedPool>，里面一把 Mutex 保护每种 dtype 一张 BTreeMap<容量, Vec<Vec<T>>>；守卫和回收器都只拿一个 Arc 克隆或弱引用，所以张量可以比借给它缓冲的那次运算活得久。总量超过上限（默认 100 MiB）时不再保留。matmul（contraction.rs 的 fresh）、超越函数、激活、softmax、线代的输出都走同一个出口，pool_recycle 程序每组运算之后看到池里留着一块 32 MiB：pool 1 buf / 32 MiB。

serial, limit    0 MiB: add 2^24 f64 (128 MiB out) x8: faults   8195  ms/call [13.24, 2.91, 2.93, 2.99, 2.93, 3.14, 3.06, 3.02]  pool BufferPoolStats { buffers: 0, capacity_bytes: 0 }
serial, limit  512 MiB: add 2^24 f64 (128 MiB out) x8: faults      0  ms/call [3.06, 2.9, 2.97, 2.94, 2.9, 3.11, 2.97, 2.9]  pool BufferPoolStats { buffers: 1, capacity_bytes: 134217728 }
serial, limit  100 MiB: add 2^24 f64 (128 MiB out) x8: faults      0  ms/call [3.16, 3.01, 2.99, 3.17, 3.04, 3.17, 3.37, 3.15]  pool BufferPoolStats { buffers: 0, capacity_bytes: 0 }

运行期在中间结果最后一次被用掉之后调 TensorBuffer::reclaim_buffer（reclaim_tensor），把它直接放回池，不等回收器；时机见执行与缓存。

融合：只在 traced 运行期发生

TensorFusion::execute_elementwise_fusion 只有一个调用方：tenferro-runtime 的 execution.rs。准备阶段 plan_elementwise_regions 把一串同 dtype 的逐元素指令划成区域，执行时一个区域换成一次融合调用。eager 和直接调用都是一个运算一个运算地执行，不会融合。

CPU engine 把区域交给 elementwise_fusion_with_pool，它返回 Ok(None) 表示「不融合」，运行期就逐条执行区域里的指令。条件依次是：

1. 输入不超过 ERASED_FUSION_MAX_INPUTS = 4 个；
2. 没有 Remainder 和 Erf（strided-fused 没有这两条指令，plan_uses_unfused_op）；
3. 不是「一个带广播的乘法」这种单指令区域，那种交给专门的广播乘法内核；
4. dtype 是 F32 或 F64（select_fused_dtype，注释说其他 dtype 只能走「the unmeasured per-element interpreter」）；
5. 所有输入（展开广播后）形状相同；
6. 元素数不少于 ELEMENTWISE_FUSION_MIN_ELEMENTS = 16 × 1024。

第 1 条在准备阶段就生效：CPU 的运行期适配器通过 max_fused_region_inputs 报告 4，规划器不会划出更大的区域（#1968，首次进入 v0.8.0）。0.7.1 里规划器照划不误，执行时报 unsupported arity 5; maximum supported arity is 4。

满足条件的区域编成一个 strided-fused 的计划，每个输出从池里取一块 PooledUninitOutput。有几个输出的区域，execute_erased_fused_plan 把指令表重放一次、写满所有输出（0.7.1 是每个输出各算一遍整个 DAG）。

strided-fused：特化、分块重放、解释器

strided-fused 0.9 按下面的顺序找执行方式（fused_elementwise_into_uninit_multi）：

1. 静态特化，只对单输出的区域：classify_static_specialization 认得单个一元、二元、三元运算，(a + b) * a 和 a * (a + b)，exp(a * b + c)，rsqrt(sqrt(clamp(a / b, lo, hi)))，命中就用单态化的 zip_map* 内核。
2. 分块寄存器重放：块里每个输出步长为 1、每个输入步长为 1 或 0（广播）时，按 1024 个元素一块（FUSED_CHUNK_ELEMS）把输入读进寄存器缓冲，每条指令对整块跑一遍单态化的内核，最后整块写出。中间结果留在 L1 里。Exp、Log、Sin、Cos、Tanh、Expm1 用的是和非融合路径同一套 elem-simd-math 内核，所以融合与不融合的结果在这几条指令上一致；chunk.rs 的注释说这条路是为了解决解释器「is slower than dispatching the same instructions one at a time」（strided-rs#304）。
3. 逐元素解释器：其余情况（f32/f64 以外的类型、步长不规则的块），每个元素过一遍寄存器数组和 match。

整个区域按 plan_class 计价：含任何超越函数就是 Transcendental，否则 Map。0.7.1 时串行上下文连静态特化也跳过、一律解释执行，现在串行令牌和并行令牌走同样的选择。

实验：一串逐元素运算

fusion 程序对四条链分别比较「在 engine 上一个个调」（direct）和「traced 编译、prepare 一次、run 反复执行」（traced），f64，串行 engine 和 8 线程 engine，再用两个 #[doc(hidden)] 的诊断方法打印区域数和「融合了几次、退回了几次」。max|diff| 是两条路结果的最大差：

(x+y)*x                n     8192 width  1: direct    0.005 ms  traced    0.007 ms  ratio  1.51  regions (1, 2) fused/fallback (0, 203)  max|diff| 0.0e0
(x+y)*x                n    16384 width  1: direct    0.009 ms  traced    0.007 ms  ratio  0.75  regions (1, 2) fused/fallback (203, 0)  max|diff| 0.0e0
(x+y)*x                n  1048576 width  1: direct    0.585 ms  traced    0.264 ms  ratio  0.45  regions (1, 2) fused/fallback (23, 0)  max|diff| 0.0e0
(x+y)*x                n  1048576 width  8: direct    0.095 ms  traced    0.074 ms  ratio  0.77  regions (1, 2) fused/fallback (23, 0)  max|diff| 0.0e0
exp(x*y+z)             n  1048576 width  1: direct    2.384 ms  traced    2.532 ms  ratio  1.06  regions (1, 3) fused/fallback (23, 0)  max|diff| 8.9e-16
exp(x*y+z)             n  1048576 width  8: direct    1.072 ms  traced    0.543 ms  ratio  0.51  regions (1, 3) fused/fallback (23, 0)  max|diff| 8.9e-16
tanh(sin(x*x+x)*y+y)   n    16384 width  1: direct    0.165 ms  traced    0.167 ms  ratio  1.01  regions (1, 6) fused/fallback (203, 0)  max|diff| 0.0e0
tanh(sin(x*x+x)*y+y)   n    16384 width  8: direct    0.166 ms  traced    0.050 ms  ratio  0.30  regions (1, 6) fused/fallback (203, 0)  max|diff| 0.0e0
tanh(sin(x*x+x)*y+y)   n  1048576 width  1: direct   11.158 ms  traced   10.577 ms  ratio  0.95  regions (1, 6) fused/fallback (23, 0)  max|diff| 0.0e0
tanh(sin(x*x+x)*y+y)   n  1048576 width  8: direct    1.854 ms  traced    1.656 ms  ratio  0.89  regions (1, 6) fused/fallback (23, 0)  max|diff| 0.0e0
a1+a2+a3+a4+a5         n  1048576 width  1: direct    1.107 ms  traced    1.128 ms  ratio  1.02  regions (0, 0) fused/fallback (0, 0)  max|diff| 0.0e0
a1+a2+a3+a4+a5         n  1048576 width  8: direct    0.516 ms  traced    0.552 ms  ratio  1.07  regions (0, 0) fused/fallback (0, 0)  max|diff| 0.0e0

（完整输出还有每条链 8192、16384 个元素的行；每格 2 次预热，2^20 时取 20 次、其余 200 次的最小值。）

• 门槛看得见：8192 个元素时三个区域都规划了，但 203 次运行全部退回逐条执行，traced 比 direct 慢 1.07 到 1.51 倍（准备好的程序每次运行的固定开销）；16384 起全部融合。
• 单线程不再吃亏：(x+y)*x 命中静态特化，2^20 个元素时融合版快 2.2 倍；exp(x*y+z) 也是特化，单线程和逐条持平；六步链走分块重放，持平或略快。同机重跑 0.7.1 的同一组实验（cpu-old 里的 fusion3、fusion2，with_threads(1)），前两条链融合后反而慢 12.5 倍和 4.7 倍（6.961 ms 对 0.556 ms，13.034 ms 对 2.791 ms）。
• 融合区域按 Transcendental 计价：六步链 16384 个元素时，8 线程的 direct 和串行一样（0.166 对 0.165 ms，五个运算都没有分派，上面问答里的原因），融合版分派了，快 3.3 倍。
• 5 个输入不再报错：规划器没有划区域（regions (0, 0)），逐条执行，结果正确。
• 8 线程、2^20 个元素的几行每轮之间波动较大（前一轮负载 10 时 (x+y)*x 的 direct 是 0.160 ms、六步链 traced 是 2.375 ms），只能看方向。
• exp 那条链有 8.9e-16 的差：静态特化 exp(a * b + c) 调的是 FusedScalar::fused_exp，即标准库的 exp，逐条执行用的是 SIMD 的 exp，两个实现差在最后一位。这也是它单线程时没有变快的原因之一：标量 exp 比 SIMD 的慢。

动脑筋答案

(x+y)*x                n  1048576 width  1: direct    0.585 ms  traced    0.264 ms  ratio  0.45  regions (1, 2) fused/fallback (23, 0)  max|diff| 0.0e0

• 融合的快，约 2.2 倍（0.264 ms 对 0.585 ms）。
• (x + y) * x 正是 classify_static_specialization 里的 AddMulLeft：一次遍历读 x、y，写一次输出，内循环单态化、可以向量化。不融合时 add 先写出一个 8 MiB 的中间结果，mul 再把它读回来，多了一整轮内存读写；单线程时这类运算受内存带宽限制，省掉一轮就是一半的时间。
• 结果逐位相同（max|diff| 0.0e0）：两条路做的是同样的两次舍入。

小结

• 一次逐元素运算在 tenferro 里按 DType 分派到泛型函数，再交给 strided-kernel 的类型擦除入口按 KernelDType 和运算码二次分派；输出永远是紧凑列主序。
• 并行由代价模型决定：「每元素耗时 × 元素数」不少于 3 倍分派代价才分派，add 在本机 131072 个元素起并行。超越函数被按 Map 计价，tanh 到 262144 个元素才并行，看起来是上游的遗漏。
• exp、tanh 等走 elem-simd-math 的 SIMD 内核，内核在 strided 的并行叶子里调用；在 M2 Max 上 f64 的 exp 比 macOS libm 快，tanh 慢将近一倍。
• 输出由 PooledUninitOutput 从池里取、不清零，assume_init_recycled 挂回收器，张量 drop 时回池，没写完的也回池。缓冲回池省掉的缺页是 glibc 上的事；macOS 的分配器自己就会复用大块内存。
• 融合只发生在 traced 运行期，要求不多于 4 个输入、没有 erf / 取余、F32 或 F64、同形状、不少于 16384 个元素；更多输入的区域不再规划，也就不再报错。
• strided-fused 0.9 先试静态特化，再试分块寄存器重放，最后才是逐元素解释器；单线程时融合不再变慢，认得的形状快约 2 倍。