📅 创建时间：2026-06-03 🏷️ 标签：#CPU #AVX #AVX512 #NEON #SIMD #向量化 #多线程 #TBB #OpenMP #Blocking #Tiling 📚 前置知识：[[13-codegen-architecture]]（代码生成架构） 📚 相关知识：[[12-scheduling]]（调度） [[17-kernel-primer]]（Kernel 开发入门）

---

CPU 后端：SIMD、分块与多线程 / CPU Backends with SIMD, Tiling, and Multithreading

┌──────────────────────────────────────────────────────────────────────────────┐ │ 情 境 描 述 │ ├──────────────────────────────────────────────────────────────────────────────┤ │ 你在 Intel Xeon Platinum 8280 (Cascade Lake) 上跑推理，发现： │ │ - 单核 1 个 inference 请求：延迟 50ms │ │ - 开 8 核并发 8 个请求：总延迟 300ms（理想应该 50ms 并行） │ │ - 你怀疑是 AVX-512 向量化没有生效，导致计算太慢，CPU 在等内存 │ └──────────────────────────────────────────────────────────────────────────────┘

第1节：CPU 执行模型——超标量流水线的奥秘

1.1 现代 CPU 架构概述

现代 CPU 是高度复杂的超标量（Superscalar）乱序执行处理器：

┌─────────────────────────────────────────────────────────────────────────────┐
│                           现代 CPU 架构示意                                  │
├─────────────────────────────────────────────────────────────────────────────┤
│                                                                             │
│                          Fetch/Decode                                       │
│                              ↓                                              │
│  ┌─────────────────────────────────────────────────────────────────────┐   │
│  │                    Out-of-Order Execution Engine                    │   │
│  │                                                                      │   │
│  │   ┌─────────┐  ┌─────────┐  ┌─────────┐  ┌─────────┐              │   │
│  │   │ ALU 0   │  │ ALU 1   │  │ FMA 0   │  │ FMA 1   │   ← 执行单元  │   │
│  │   └────┬────┘  └────┬────┘  └────┬────┘  └────┬────┘              │   │
│  │        └──────┬─────┴──────┬─────┴──────┬─────┘                    │   │
│  │               ↓             ↓            ↓                            │   │
│  │         ┌─────────────────────────────────────────┐                  │   │
│  │         │         Register File (32/64 GPR)       │                  │   │
│  │         └──────────────────┬──────────────────────┘                  │   │
│  │                            ↓                                          │   │
│  │         ┌─────────────────────────────────────────┐                  │   │
│  │         │              Reorder Buffer              │                  │   │
│  │         └──────────────────┬──────────────────────┘                  │   │
│  └────────────────────────────┼──────────────────────────────────────────┘   │
│                               ↓                                              │
│                          Retire (提交)                                        │
│                                                                             │
└─────────────────────────────────────────────────────────────────────────────┘

1.2 CPU 流水线的关键概念

class CPUPipeline:
    """
    CPU 流水线关键概念
    """
    
    # 1. 流水线级数 (Pipeline Depth)
    # Intel Skylake: 14 级流水线
    pipeline_stages = {
        'Skylake': 14,      # Fetch→Decode→Rename→Dispatch→Execute→Retire
        'Zen 2': 12,        # AMD
        'Apple M1': 12      # ARM custom
    }
    
    # 2. 每周期发射数 (IPC - Instructions Per Cycle)
    # Skylake: 4-wide issue (每周期最多 4 条指令发射到执行单元)
    max_issue_width = 4
    
    # 3. 执行单元数量
    execution_units = {
        'ALU': 4,           # 算术逻辑单元
        'FMA': 2,           # 乘加单元（用于 FP 计算）
        'Load': 2,          # 内存加载
        'Store': 1          # 内存存储
    }
    
    # 4. 分支预测
    # 采用 TAGE (Tagged Geometric History Length) 预测器
    # 准确率 > 95%
    branch_predictor = 'TAGE-SC-L'
    
    # 5. 乱序执行窗口 (ROB - Reorder Buffer)
    # Skylake: 224 项
    rob_size = 224
    
    # 6. 寄存器重命名
    # Architectural GPR: 16 (RAX, RBX, ...)
    # Physical Registers: 180 (Skylake)
    arch_regs = 16
    phys_regs = 180


def pipeline_stall_analysis():
    """
    流水线停顿分析
    """
    
    # 常见停顿原因：
    
    # 1. 数据依赖停顿 (Data Hazard)
    # 必须等待前一条指令的结果
    code = """
        vfmadd231ps zmm0, zmm1, zmm2   ; zmm0 = zmm0 + zmm1 * zmm2
        vfmadd231ps zmm0, zmm3, zmm4   ; 必须等待上面的 zmm0
        vfmadd231ps zmm0, zmm5, zmm6   ; 必须等待上面的 zmm0
    """
    # 依赖链长 = 3 条指令 × 4 cycles (FMA 延迟) = 12 cycles stall
    
    # 2. 分支预测失败 (Branch Misprediction)
    # Skylake: 约 15 cycles 惩罚
    mispredict_penalty = 15
    
    # 3. 缓存未命中 (Cache Miss)
    # L1: ~4 cycles
    # L2: ~12 cycles
    # L3: ~40 cycles
    # DRAM: ~200 cycles
    cache_latency = {'L1': 4, 'L2': 12, 'L3': 40, 'DRAM': 200}
    
    # 4. 资源冲突
    # 所有执行单元被占用
    resource_stall = 0  # 4-wide issue 很少成为瓶颈

1.3 CPU 性能公式

def cpu_performance_formula():
    """
    CPU 性能三大要素
    """
    
    # 公式：执行时间 = 指令数 × CPI × 时钟周期
    
    # 1. 指令数 (IC - Instruction Count)
    # 取决于算法和编译器
    # 向量化可以将多条标量指令合并为一条向量指令
    ic_scalar = 1000000  # 标量版本
    ic_vector = 125000    # 向量化版本（假设 SIMD 宽度 = 8）
    
    # 2. 每指令周期数 (CPI - Cycles Per Instruction)
    # 理想情况：1 (每周期完成一条指令)
    # 实际受 stall 和依赖影响
    cpi_ideal = 1.0
    cpi_with_stalls = 2.5  # 考虑内存延迟
    
    # 3. 时钟周期
    # Intel Xeon 8280: 2.5 GHz
    clock_ghz = 2.5
    
    # 计算
    time_scalar = ic_scalar * cpi_with_stalls / (clock_ghz * 1e9)
    time_vector = ic_vector * cpi_with_stalls / (clock_ghz * 1e9)
    
    speedup = time_scalar / time_vector
    print(f"向量化加速比: {speedup:.1f}x")

第2节：SIMD 基础——一次计算多个数据

2.1 SIMD 指令集演进

┌─────────────────────────────────────────────────────────────────────────────┐
│                           SIMD 指令集演进                                    │
├─────────────────────────────────────────────────────────────────────────────┤
│                                                                             │
│  x86 (Intel/AMD):                                                          │
│  ┌──────────┬───────────┬───────────┬─────────────┬───────────────┐       │
│  │  SSE     │  SSE2/3/4 │  AVX      │  AVX2       │  AVX-512      │       │
│  │  (1999)  │  (2000s)  │  (2011)   │  (2013)     │  (2016)       │       │
│  ├──────────┼───────────┼───────────┼─────────────┼───────────────┤       │
│  │ 128-bit  │  128-bit  │  256-bit  │  256-bit    │  512-bit      │       │
│  │ 4× FP32  │  + int    │  8× FP32  │  32 byte ops│  16× FP32     │       │
│  │          │           │  4× FP64  │  + gather   │  8× FP64      │       │
│  │          │           │           │             │  32× int8     │       │
│  └──────────┴───────────┴───────────┴─────────────┴───────────────┘       │
│                                                                             │
│  ARM:                                                                      │
│  ┌──────────┬───────────┬───────────┬─────────────┬───────────────┐       │
│  │  NEON    │  SVE      │  SVE2     │  SME        │               │       │
│  │  (2005)  │  (2016)   │  (2019)   │  (2022)     │               │       │
│  ├──────────┼───────────┼───────────┼─────────────┼───────────────┤       │
│  │  128-bit │  可变长度  │  可变长度 │  可变长度   │               │       │
│  │  SIMD    │  2-2048   │  +更多    │  可扩展     │               │       │
│  │          │  元素      │  指令     │  矩阵运算   │               │       │
│  └──────────┴───────────┴───────────┴─────────────┴───────────────┘       │
│                                                                             │
└─────────────────────────────────────────────────────────────────────────────┘

2.2 SIMD 宽度与数据格式

指令集	宽度	FP32 元素	FP64 元素	INT8 元素	主要应用
SSE	128-bit	4	2	16	遗留代码
AVX	256-bit	8	4	32	通用计算
AVX2	256-bit	8	4	32	通用计算 + Gather
AVX-512	512-bit	16	8	64	高性能计算
AVX-512 VNNI	512-bit	-	-	64	INT8 推理
AVX-512 BF16	512-bit	-	8*	64	BFloat16 推理
NEON	128-bit	4	2	16	ARM 移动端
SVE	可变	N	N/2	2N	ARM HPC

2.3 Intrinsics 编程基础

// AVX-512 Intrinsics 示例：向量加法
#include <immintrin.h>

void vector_add_avx512(float* c, const float* a, const float* b, int n) {
    // c[i] = a[i] + b[i], 向量化版本
    for (int i = 0; i < n; i += 16) {  // 16 = 512bit / 32bit
        // 加载 16 个 float (512 bits)
        __m512 va = _mm512_loadu_ps(&a[i]);
        __m512 vb = _mm512_loadu_ps(&b[i]);
        
        // 向量加法
        __m512 vc = _mm512_add_ps(va, vb);
        
        // 存储结果
        _mm512_storeu_ps(&c[i], vc);
    }
    
    // 处理尾部（n 可能不是 16 的倍数）
    for (int i = (n / 16) * 16; i < n; i++) {
        c[i] = a[i] + b[i];
    }
}

// AVX-512 FMA：融合乘加
void vector_gemm_element_avx512(float* c, const float* a, const float* b, int K) {
    // 计算 c += sum_k(A[i,k] * B[k,j])
    // 使用 FMA 指令：d = a * b + c
    __m512 vc = _mm512_setzero_ps();  // 初始化为 0
    
    for (int k = 0; k < K; k++) {
        __m512 va = _mm512_set1_ps(a[k]);     // 广播 a[k]
        __m512 vb = _mm512_loadu_ps(&b[k * 16]);  // 加载 b[k, 0:16]
        vc = _mm512_fmadd_ps(va, vb, vc);    // vc = vc + va * vb
    }
    
    _mm512_storeu_ps(c, vc);
}

// AVX-512 Masked 操作：处理动态长度
void vector_add_masked_avx512(float* c, const float* a, const float* b, int n) {
    for (int i = 0; i < n; i += 16) {
        int remaining = n - i;
        
        if (remaining >= 16) {
            // 完整向量
            __m512 va = _mm512_loadu_ps(&a[i]);
            __m512 vb = _mm512_loadu_ps(&b[i]);
            __m512 vc = _mm512_add_ps(va, vb);
            _mm512_storeu_ps(&c[i], vc);
        } else {
            // 部分向量，使用 mask
            // 创建 mask: 低 remaining 位为 1
            __mmask16 mask = (1 << remaining) - 1;
            
            __m512 va = _mm512_maskz_loadu_ps(mask, &a[i]);
            __m512 vb = _mm512_maskz_loadu_ps(mask, &b[i]);
            __m512 vc = _mm512_add_ps(va, vb);
            _mm512_mask_storeu_ps(&c[i], mask, vc);
        }
    }
}

2.4 Python NumPy 对比

import numpy as np

def numpy_vs_intrinsics():
    """
    NumPy (使用 MKL/OpenBLAS) vs 手写 Intrinsics
    """
    
    n = 1000000
    a = np.random.randn(n).astype(np.float32)
    b = np.random.randn(n).astype(np.float32)
    
    # NumPy 版本（会自动使用 SIMD）
    start = time.time()
    c_numpy = a + b
    time_numpy = time.time() - start
    
    # 手写 Intrinsics 需要 C 扩展
    # 这里假设已经编译成 .so
    c_c = vector_add_avx512_c(a, b)  # 伪代码
    
    print(f"NumPy: {time_numpy*1000:.2f} ms")
    print(f"手写 Intrinsics: {time_c*1000:.2f} ms")
    
    # 理论上性能相近，因为 MKL 也是用 Intrinsics 写的


def numpy_gemm():
    """
    NumPy GEMM 性能分析
    """
    
    M, K, N = 1024, 512, 1024
    A = np.random.randn(M, K).astype(np.float32)
    B = np.random.randn(K, N).astype(np.float32)
    
    # NumPy matmul 使用 MKL/OpenBLAS，会自动：
    # 1. 选择合适的 blocking 大小
    # 2. 使用 SIMD 向量化
    # 3. 使用多线程
    start = time.time()
    C = A @ B
    time_gemm = time.time() - start
    
    gflops = 2 * M * K * N / time_gemm / 1e9
    print(f"GEMM: {time_gemm*1000:.2f} ms, {gflops:.1f} GFLOPS")

第3节：向量化——自动 vs 手动

3.1 自动向量化

现代编译器可以自动向量化简单循环：

// 原始代码
void add(float* c, const float* a, const float* b, int n) {
    for (int i = 0; i < n; i++) {
        c[i] = a[i] + b[i];
    }
}

// GCC -O3 自动向量化（检查：gcc -O3 -ftree-vectorize -S -fopt-info-vec-optimized）
// 编译器会生成 SIMD 指令
// 使用 -fopt-info-vec-optimized 查看向量化信息

编译器自动向量化的条件：

条件	说明
循环次数可确定	必须是固定上界或可推断
数据对齐	指针 32/64 字节对齐更好
无依赖	循环迭代之间无数据依赖
简单操作	编译器能识别的模式
启用优化	-O2 或 -O3

自动向量化的局限性：

// 无法自动向量化的例子
void problematic(float* a, float* b, int n) {
    for (int i = 1; i < n; i++) {
        a[i] = a[i-1] + b[i];  // 循环携带依赖，无法向量化
    }
}

// 无法识别的模式
void complex_pattern(float* a, int n) {
    for (int i = 0; i < n; i++) {
        if (a[i] > 0) {  // 条件分支
            a[i] = sqrt(a[i]);
        } else {
            a[i] = -sqrt(-a[i]);
        }
    }
}

3.2 手动向量化

当自动向量化不够时，需要手动控制：

// GEMM 手动向量化
// c[i,j] += sum_k(A[i,k] * B[k,j])
// 使用 blocking + 手动 SIMD

void sgemm_blocked_avx512(
    float* C, const float* A, const float* B,
    int M, int N, int K,
    int block_m, int block_n, int block_k) {
    
    // 外层 blocking
    for (int jj = 0; jj < N; jj += block_n) {
        for (int kk = 0; kk < K; kk += block_k) {
            // 内层 blocking
            for (int i = 0; i < M; i += block_m) {
                
                // 寄存器级别的 K 循环
                int j_end = min(jj + block_n, N);
                int k_end = min(kk + block_k, K);
                int i_end = min(i + block_m, M);
                
                // 累加到 C[i:i_end, jj:j_end]
                for (int k = kk; k < k_end; k++) {
                    // 加载 A[i:i_end, k]（一列）
                    // 使用广播
                    for (int ii = i; ii < i_end; ii += 16) {
                        __m512 a_vec = _mm512_loadu_ps(&A[ii * K + k]);
                        
                        // 对 B[k, jj:j_end] 的每个 16-element 块
                        for (int j = jj; j < j_end; j += 16) {
                            __m512 b_vec = _mm512_loadu_ps(&B[k * N + j]);
                            __m512 c_vec = _mm512_loadu_ps(&C[ii * N + j]);
                            c_vec = _mm512_fmadd_ps(a_vec, b_vec, c_vec);
                            _mm512_storeu_ps(&C[ii * N + j], c_vec);
                        }
                    }
                }
            }
        }
    }
}

3.3 性能对比

方法	FP32 峰值	相对性能	适用场景
纯 Python	-	1x	原型开发
NumPy (单核)	~50 GFLOPS	~10x	一般使用
MKL 单核	~200 GFLOPS	~40x	需要极致单核
MKL 多核	~1600 GFLOPS	~320x	生产部署
手写 Intrinsics	~150 GFLOPS	~30x	学习/特殊优化
手写 + 多核	~1200 GFLOPS	~240x	最优性能

第4节：内存访问优化——Cache 为王

4.1 CPU Cache 层级

┌─────────────────────────────────────────────────────────────────────────────┐
│                           CPU Cache 层级                                    │
├─────────────────────────────────────────────────────────────────────────────┤
│                                                                             │
│                              CPU Core                                       │
│  ┌───────────────────────────────────────────────────────────────────────┐ │
│  │  Registers (R0-R31):  32 × 512-bit = 2 KB         ← 最快，最贵       │ │
│  │  Latency: 1 cycle                                              ↑      │ │
│  └───────────────────────────────────────────────────────────────────────┘ │
│                                      ↓ Latency                              │
│  ┌───────────────────────────────────────────────────────────────────────┐ │
│  │  L1 Instruction Cache: 32 KB                                     ↑      │ │
│  │  L1 Data Cache: 32 KB                                         ↑        │ │
│  │  Latency: 4 cycles                                            ↑          │ │
│  └───────────────────────────────────────────────────────────────────────┘ │
│                                      ↓ Latency                              │
│  ┌───────────────────────────────────────────────────────────────────────┐ │
│  │  L2 Unified Cache: 256 KB - 1 MB                              ↑        │ │
│  │  Latency: 12 cycles                                           ↑          │ │
│  └───────────────────────────────────────────────────────────────────────┘ │
│                                      ↓ Latency                              │
│  ┌───────────────────────────────────────────────────────────────────────┐ │
│  │  L3 Shared Cache: 8 MB - 64 MB (per socket)                    ↑        │ │
│  │  Latency: 40 cycles                                           ↑          │ │
│  └───────────────────────────────────────────────────────────────────────┘ │
│                                      ↓ Latency                              │
│  ┌───────────────────────────────────────────────────────────────────────┐ │
│  │  Main Memory (DDR4/DDR5): TB                                     ↑        │ │
│  │  Latency: 200 cycles                                           ↑          │ │
│  └───────────────────────────────────────────────────────────────────────┘ │
│                                                                             │
└─────────────────────────────────────────────────────────────────────────────┘

4.2 Blocking / Tiling

Blocking 是最重要的内存优化技术：

// 未优化的 GEMM
// 计算: C = A × B
// 复杂度: O(M × N × K)
// 内存访问: C 的每个元素读取 K 次
void gemm_naive(float* C, const float* A, const float* B, int M, int N, int K) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
    // 问题: A[i,k] 被读取 M 次！
    //       B[k,j] 被读取 N 次！
    //       Cache 效率极低
}

// Blocked GEMM - 将大矩阵分块
// 核心思想: 将 A, B, C 分成小块，每个小块可以放入 L1 Cache
void gemm_blocked(float* C, const float* A, const float* B, 
                  int M, int N, int K, int block_size) {
    
    // 分配 L1 Cache 大小的块
    // 假设 block_size = 64, 则每个块 64×64×4 = 16 KB < L1 Cache 32 KB
    
    for (int i0 = 0; i0 < M; i0 += block_size) {
        for (int j0 = 0; j0 < N; j0 += block_size) {
            // 初始化 C 块为 0
            for (int i = i0; i < min(i0 + block_size, M); i++) {
                for (int j = j0; j < min(j0 + block_size, N); j++) {
                    C[i * N + j] = 0.0f;
                }
            }
            
            // K 维度分块
            for (int k0 = 0; k0 < K; k0 += block_size) {
                // 加载 A 块到 L1 Cache
                // A_block[i, k] = A[i, k0 + k]
                
                // 加载 B 块到 L1 Cache
                // B_block[k, j] = B[k0 + k, j]
                
                // 计算 C 块 += A 块 × B 块
                for (int i = i0; i < min(i0 + block_size, M); i++) {
                    for (int k = k0; k < min(k0 + block_size, K); k++) {
                        float a_ik = A[i * K + k];
                        
                        // 内层 j 循环
                        // 这个循环可以被编译器向量化
                        for (int j = j0; j < min(j0 + block_size, N); j++) {
                            C[i * N + j] += a_ik * B[k * N + j];
                        }
                    }
                }
            }
        }
    }
}

4.3 Loop Reordering

改变循环顺序可以改善局部性：

// 原始顺序（列优先访问 A，行优先访问 B）
// A 按列访问，空间局部性差
for (int k = 0; k < K; k++) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            C[i * N + j] += A[i * K + k] * B[k * N + j];
        }
    }
}

// 重排后的顺序（行优先访问 A，列优先访问 B）
// A 和 B 都按行访问，Cache 友好
for (int i = 0; i < M; i++) {
    for (int k = 0; k < K; k++) {
        float a_ik = A[i * K + k];
        for (int j = 0; j < N; j++) {
            C[i * N + j] += a_ik * B[k * N + j];
        }
    }
}

// 最佳顺序：i-j-k 循环顺序，k 作为最内层循环
// 使得 A[i,k] 和 B[k,j] 都被连续访问
for (int i = 0; i < M; i++) {
    for (int j = 0; j < N; j++) {
        float sum = 0.0f;
        for (int k = 0; k < K; k++) {
            sum += A[i * K + k] * B[k * N + j];
        }
        C[i * N + j] = sum;
    }
}

4.4 Prefetch

Prefetch 提前将数据加载到 Cache：

// 手动 Prefetch
void gemm_with_prefetch(float* C, const float* A, const float* B, 
                        int M, int N, int K, int block_size) {
    
    for (int i0 = 0; i0 < M; i0 += block_size) {
        for (int j0 = 0; j0 < N; j0 += block_size) {
            // 初始化 C 块
            for (int i = i0; i < min(i0 + block_size, M); i++) {
                for (int j = j0; j < min(j0 + block_size, N); j++) {
                    C[i * N + j] = 0.0f;
                }
            }
            
            for (int k0 = 0; k0 < K; k0 += block_size) {
                // 加载 A 块和 B 块
                
                // Prefetch 下一块 B（当 k 很大时，数据可能已经换出）
                if (k0 + block_size < K) {
                    for (int j = j0; j < min(j0 + block_size, N); j++) {
                        // _mm_prefetch(&B[(k0 + block_size) * N + j], _MM_HINT_T0);
                        // 提示 CPU 将数据预取到 L1 Cache
                    }
                }
                
                // 计算
                for (int i = i0; i < min(i0 + block_size, M); i++) {
                    for (int k = k0; k < min(k0 + block_size, K); k++) {
                        float a_ik = A[i * K + k];
                        for (int j = j0; j < min(j0 + block_size, N); j++) {
                            C[i * N + j] += a_ik * B[k * N + j];
                        }
                    }
                }
            }
        }
    }
}

第5节：多线程并行——TBB vs OpenMP

5.1 OpenMP 编程

// OpenMP GEMM
#include <omp.h>

void gemm_omp(float* C, const float* A, const float* B, 
              int M, int N, int K, int num_threads) {
    
    omp_set_num_threads(num_threads);
    
    #pragma omp parallel for collapse(2) schedule(static)
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            float sum = 0.0f;
            for (int k = 0; k < K; k++) {
                sum += A[i * K + k] * B[k * N + j];
            }
            C[i * N + j] = sum;
        }
    }
}

// 编译：gcc -O3 -fopenmp gemm.c -o gemm
// 设置线程数：export OMP_NUM_THREADS=8

5.2 Intel TBB 编程

// Intel TBB GEMM
#include <tbb/tbb.h>
#include <tbb/blocked_range2d.h>

using namespace tbb;

class GemmBody {
public:
    float* C;
    const float* A;
    const float* B;
    int N;
    int K;
    
    GemmBody(float* c, const float* a, const float* b, int n, int k)
        : C(c), A(a), B(b), N(n), K(k) {}
    
    void operator()(const blocked_range2d<int>& r) const {
        for (int i = r.rows().begin(); i < r.rows().end(); i++) {
            for (int j = r.cols().begin(); j < r.cols().end(); j++) {
                float sum = 0.0f;
                for (int k = 0; k < K; k++) {
                    sum += A[i * K + k] * B[k * N + j];
                }
                C[i * N + j] = sum;
            }
        }
    }
};

void gemm_tbb(float* C, const float* A, const float* B, 
              int M, int N, int K) {
    
    blocked_range2d<int> range(0, M, 0, N);
    
    parallel_for(range, GemmBody(C, A, B, N, K));
}

5.3 多线程策略对比

特性	OpenMP	Intel TBB	std::thread
语法	编译指示 (#pragma)	库函数	手动管理
负载均衡	静态/动态调度	自动	手动
任务依赖	有限	丰富	无
跨平台	好	好	好
学习曲线	低	中	高
适用场景	简单并行	复杂任务图	细粒度控制

第6节：Intel AMX——下一代 CPU 加速

6.1 AMX 架构

Intel AMX (Advanced Matrix Extensions) 是专门为矩阵计算设计的指令集：

┌─────────────────────────────────────────────────────────────────────────────┐
│                           AMX 架构                                          │
├─────────────────────────────────────────────────────────────────────────────┤
│                                                                             │
│  核心组件:                                                                   │
│  ┌───────────────────────────────────────────────────────────────────────┐   │
│  │  Tile Register File: 8 × 1024-byte tiles = 8 KB                       │   │
│  │  每个 tile 是一个 16×16 FP32/BF16 矩阵                                │   │
│  └───────────────────────────────────────────────────────────────────────┘   │
│                              ↓                                              │
│  核心指令:                                                                   │
│  ┌───────────────────────────────────────────────────────────────────────┐   │
│  │  LD_TILECFG: 配置 tile 形状                                           │   │
│  │  ST_TILECFG: 保存 tile 配置                                           │   │
│  │  LDTILEF: 加载数据到 tile (from memory)                               │   │
│  │  STTILEF: 保存 tile 数据 (to memory)                                  │   │
│  │  TDPBFSS: BF16 矩阵乘加 (Dot Product BF16)                           │   │
│  │  TDPBFSN: BF16 矩阵乘加 (with saturation)                            │   │
│  │  TDPFP32: FP32 矩阵乘加 (高精度)                                      │   │
│  └───────────────────────────────────────────────────────────────────────┘   │
│                                                                             │
│  性能:                                                                      │
│  ┌───────────────────────────────────────────────────────────────────────┐   │
│  │  每个 tile: 16×16 = 256 elements                                      │   │
│  │  每周期: 1 TDPBFSS = 256 BF16 × 2 (乘加) = 512 FLOPS/cycle             │   │
│  │  2.0 GHz × 512 FLOPS × 4 cores = 4096 GFLOPS (BF16)                   │   │
│  │  相比 AVX-512: ~16x 提升                                             │   │
│  └───────────────────────────────────────────────────────────────────────┘   │
│                                                                             │
└─────────────────────────────────────────────────────────────────────────────┘

6.2 AMX GEMM 示例

// AMX GEMM 伪代码（实际需要 intrinsics 或编译器支持）
void amx_gemm(float* C, const float* A, const float* B, 
             int M, int N, int K) {
    
    // 1. 配置 tile 形状
    // 每 tile: 16 行 (M), 16 列 (N/K)
    tile_config tc = {
        .rows = 16,
        .cols = 16,
        .shape = TILE_16x16
    };
    _tile_config(&tc);
    
    // 2. 主循环
    for (int m = 0; m < M; m += 16) {
        for (int n = 0; n < N; n += 16) {
            // 清零 C tile
            _tile_zero(0);
            
            // K 维度累加
            for (int k = 0; k < K; k += 16) {
                // 加载 A tile: C[m:m+16, k:k+16]
                _ldtilef(&A[m * K + k], 1);  // tmm0
                
                // 加载 B tile: B[k:k+16, n:n+16]
                _ldtilef(&B[k * N + n], 2);   // tmm1
                
                // C += A × B
                // tdpbfss: BF16 输入, FP32 累加
                _tdpbfss(0, 1, 2);  // dest=tmm0, src1=tmm1, src2=tmm2
            }
            
            // 保存 C tile
            _sttilef(&C[m * N + n], 0);
        }
    }
}

---

升华

┌─────────────────────────────────────────────────────────────────────────────────┐ │ CPU 优化的核心原则 │ ├─────────────────────────────────────────────────────────────────────────────────┤ │ │ │ 1. Cache 为王：所有优化的核心是让数据在 L1 Cache 中，最大化重用 │ │ │ │ 2. 向量化是基础：SIMD 指令是 CPU 的原生能力，不用白不用 │ │ │ │ 3. 多线程要谨慎：同步开销和 Cache 一致性会影响并行效率 │ │ │ │ 4. AMX 是未来：BF16 + Tile 指令将 CPU 矩阵计算提升到新水平 │ │ │ └─────────────────────────────────────────────────────────────────────────────────┘

"AI 可查 vs 必须理解"清单

必须理解（不理解就等于不会）：

• 🔴 L1/L2/L3 Cache 的延迟差异——不知道就无法理解为什么要 blocking
• 🔴 SIMD 向量化的原理——一条指令操作多个数据，不知道就无法理解性能来源
• 🔴 Blocking / Tiling 的动机——将大矩阵分块适应 Cache 容量
• 🔴 手动向量化与自动向量化的边界——哪些编译器能做，哪些必须手写

AI 可查（知道去哪查就行）：

• ✅ 具体 CPU 的 Cache 大小（如 Xeon 8280 的 L3 = 38.5 MB）——查 Intel ARK
• ✅ AVX-512 Intrinsics 的具体 API（如 _mm512_fmadd_ps 的参数）——查 Intel Intrinsics Guide
• ✅ AMX 指令的具体语法——需要 Intel 官方文档
• ✅ OpenMP/TBB 的最新 API 变化——库版本相关

学习状态：🟡 开始学习