Single Instruction, Multiple Data (SIMD)

本文为学习 Computer Enhance - Single Instruction, Multiple Data 的整理。

回顾与引入

在前一篇博客中，我们通过 打破指令依赖链，让 CPU 可以并行执行 add 操作，从而提升了性能。

在这篇博客中，我们回到第一种优化方式：通过使用 SIMD 指令来减少指令总数。

SIMD 全称 Single Instruction, Multiple Data（单指令多数据），它的意思正如字面所示：CPU 用一条指令同时处理多份数据。

---

SIMD 基础

常见的 SIMD 指令集包括 SSE、AVX 和 AVX-512。指令宽度依次加倍，分别是 128 bit、256 bit 和 512 bit。宽度决定了一条指令一次性可以计算的数据长度。

现在的家用电脑一般都支持到 AVX2。下图是我的 CPU 数据，可以看到对 SIMD 指令集的支持情况：

寄存器对应关系

SIMD 指令会使用专用的寄存器做计算：

• SSE (128位): 使用 XMM0 ~ XMM15
• AVX (256位): 使用 YMM0 ~ YMM15
  • XMM0 是 YMM0 的低 128 位
• AVX-512 (512位): 使用 ZMM0 ~ ZMM31
  • YMM0 是 ZMM0 的低 256 位

---

汇编实例：8 个 float 加法

假设有两个数组 A 和 B 在内存里，想把结果存到 C。

步骤 A：搬运 (Load)

先把数据从内存搬到寄存器。

vmovups ymm0, [地址A]  ; 把 A 的 8 个 float 搬进 YMM0
vmovups ymm1, [地址B]  ; 把 B 的 8 个 float 搬进 YMM1

步骤 B：计算 (Calculate)

CPU 内部的加法器开始工作。

vaddps ymm2, ymm0, ymm1 ; 结果 = YMM0 + YMM1，放入 YMM2

• v = Vector (AVX前缀)
• add = 加法
• ps = Packed Single (打包的单精度浮点)

步骤 C：存回 (Store)

把结果搬回内存。

vmovups [地址C], ymm2  ; 把 YMM2 的结果写回内存 C

---

为什么 SIMD 更快？

使用一条 SIMD 指令 vs 使用几条普通指令，区别在哪里？为什么可以达到巨大的优化效果？

SIMD 并不能减少 CPU ALU（计算单元）的工作量，但它极大地减少了前端的取指译码以及后端的调度与寄存器访问开销。

与其让 CPU 为 4 个独立的加法做 4 套全流程管理，不如用 1 条指令打包处理，让管理成本降到最低。

具体开销对比如下：

1. 前端 (Front-end)

• 取指 (Fetch)
  • 标量：占用 4 个取指名额。瞬间吃光每周期的取指带宽（通常仅 4-6 条），导致后续指令进不来。
  • SIMD：只占用 1 个 取指名额。节省 75% 的带宽。
• 译码 (Decode)
  • 标量：进行 4 次译码。译码器是最复杂、耗电的单元。
  • SIMD：只进行 1 次 译码。

2. 后端 (Back-end)

• ROB (重排序缓冲区)
  • 标量：占 4 个格子。ROB 满了 CPU 就会阻塞。
  • SIMD：只占 1 个 格子。节省 75% 空间。
• 发射队列 (Issue Queue)
  • 标量：每周期要挑出 4 个指令，仲裁逻辑复杂。
  • SIMD：只要挑 1 个。
• 寄存器读写 (Register Ports)
  • 标量：读 8 次、写 4 次。容易卡死后端物理瓶颈。
  • SIMD：只读 2 次宽寄存器、写 1 次。
• 依赖检查
  • 标量：检查 4 条指令的依赖。
  • SIMD：只检查 1 次。

注意：使用 SIMD 并不能自动解决串行依赖链的问题，我们仍然可以通过 使用临时变量（多累加器）来打破依赖链。

---

代码演示

以下代码包含了前两章的初始代码、循环展开优化，以及本章的 SIMD 优化代码。

/* ========================================================================
   IPC (Instructions Per Clock) & SIMD Demo
   展示循环展开、打破依赖链以及 SIMD (SSE/AVX) 对性能的影响
   ======================================================================== */

#include <immintrin.h> // SIMD intrinsics
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

typedef unsigned int u32;

// ========================================================================
// Scalar Versions (IPC Focus)
// ========================================================================

// 原始的单标量累加
u32 SingleScalar(u32 Count, u32 *Input) {
  u32 Sum = 0;
  for (u32 Index = 0; Index < Count; ++Index) {
    Sum += Input[Index];
  }
  return Sum;
}

// 循环展开2次
u32 Unroll2Scalar(u32 Count, u32 *Input) {
  u32 Sum = 0;
  for (u32 Index = 0; Index < Count; Index += 2) {
    Sum += Input[Index];
    Sum += Input[Index + 1];
  }
  return Sum;
}

// 循环展开4次
u32 Unroll4Scalar(u32 Count, u32 *Input) {
  u32 Sum = 0;
  for (u32 Index = 0; Index < Count; Index += 4) {
    Sum += Input[Index];
    Sum += Input[Index + 1];
    Sum += Input[Index + 2];
    Sum += Input[Index + 3];
  }
  return Sum;
}

// 双累加器 (Breaking Dependency)
u32 DualScalar(u32 Count, u32 *Input) {
  u32 SumA = 0;
  u32 SumB = 0;
  for (u32 Index = 0; Index < Count; Index += 2) {
    SumA += Input[Index + 0];
    SumB += Input[Index + 1];
  }
  u32 Sum = SumA + SumB;
  return Sum;
}

// 四累加器
u32 QuadScalar(u32 Count, u32 *Input) {
  u32 SumA = 0;
  u32 SumB = 0;
  u32 SumC = 0;
  u32 SumD = 0;
  for (u32 Index = 0; Index < Count; Index += 4) {
    SumA += Input[Index + 0];
    SumB += Input[Index + 1];
    SumC += Input[Index + 2];
    SumD += Input[Index + 3];
  }
  u32 Sum = SumA + SumB + SumC + SumD;
  return Sum;
}

// 四累加器 + 指针递增
u32 QuadScalarPtr(u32 Count, u32 *Input) {
  u32 SumA = 0;
  u32 SumB = 0;
  u32 SumC = 0;
  u32 SumD = 0;

  Count /= 4;
  while (Count--) {
    SumA += Input[0];
    SumB += Input[1];
    SumC += Input[2];
    SumD += Input[3];
    Input += 4;
  }

  u32 Sum = SumA + SumB + SumC + SumD;
  return Sum;
}

// ========================================================================
// SIMD Versions (SSE & AVX)
// ========================================================================

// SSE Version (4-wide)
u32 __attribute__((target("ssse3"))) SingleSSE(u32 Count, u32 *Input) {
  __m128i Sum = _mm_setzero_si128();
  for (u32 Index = 0; Index < Count; Index += 4) {
    Sum = _mm_add_epi32(Sum, _mm_loadu_si128((__m128i *)&Input[Index]));
  }
  // 水平相加 reduce
  Sum = _mm_hadd_epi32(Sum, Sum);
  Sum = _mm_hadd_epi32(Sum, Sum);
  return _mm_cvtsi128_si32(Sum);
}

// AVX Version (8-wide)
u32 __attribute__((target("avx2"))) SingleAVX(u32 Count, u32 *Input) {
  __m256i Sum = _mm256_setzero_si256();
  for (u32 Index = 0; Index < Count; Index += 8) {
    Sum = _mm256_add_epi32(Sum, _mm256_loadu_si256((__m256i *)&Input[Index]));
  }
  // Reduce
  Sum = _mm256_hadd_epi32(Sum, Sum);
  Sum = _mm256_hadd_epi32(Sum, Sum);
  __m256i SumS = _mm256_permute2x128_si256(Sum, Sum, 1 | (1 << 4));
  Sum = _mm256_add_epi32(Sum, SumS);

  return _mm256_cvtsi256_si32(Sum);
}

// AVX + Unrolled 2x (Dual Accumulators)
u32 __attribute__((target("avx2"))) DualAVX(u32 Count, u32 *Input) {
  __m256i SumA = _mm256_setzero_si256();
  __m256i SumB = _mm256_setzero_si256();
  for (u32 Index = 0; Index < Count; Index += 16) {
    SumA = _mm256_add_epi32(SumA, _mm256_loadu_si256((__m256i *)&Input[Index]));
    SumB = _mm256_add_epi32(SumB, _mm256_loadu_si256((__m256i *)&Input[Index + 8]));
  }

  __m256i Sum = _mm256_add_epi32(SumA, SumB);
  // Reduce...
  Sum = _mm256_hadd_epi32(Sum, Sum);
  Sum = _mm256_hadd_epi32(Sum, Sum);
  __m256i SumS = _mm256_permute2x128_si256(Sum, Sum, 1 | (1 << 4));
  Sum = _mm256_add_epi32(Sum, SumS);

  return _mm256_cvtsi256_si32(Sum);
}

// AVX + Unrolled 4x (Quad Accumulators)
u32 __attribute__((target("avx2"))) QuadAVX(u32 Count, u32 *Input) {
  __m256i SumA = _mm256_setzero_si256();
  __m256i SumB = _mm256_setzero_si256();
  __m256i SumC = _mm256_setzero_si256();
  __m256i SumD = _mm256_setzero_si256();
  for (u32 Index = 0; Index < Count; Index += 32) {
    SumA = _mm256_add_epi32(SumA, _mm256_loadu_si256((__m256i *)&Input[Index]));
    SumB = _mm256_add_epi32(SumB, _mm256_loadu_si256((__m256i *)&Input[Index + 8]));
    SumC = _mm256_add_epi32(SumC, _mm256_loadu_si256((__m256i *)&Input[Index + 16]));
    SumD = _mm256_add_epi32(SumD, _mm256_loadu_si256((__m256i *)&Input[Index + 24]));
  }

  __m256i SumAB = _mm256_add_epi32(SumA, SumB);
  __m256i SumCD = _mm256_add_epi32(SumC, SumD);
  __m256i Sum = _mm256_add_epi32(SumAB, SumCD);
  
  // Reduce...
  Sum = _mm256_hadd_epi32(Sum, Sum);
  Sum = _mm256_hadd_epi32(Sum, Sum);
  __m256i SumS = _mm256_permute2x128_si256(Sum, Sum, 1 | (1 << 4));
  Sum = _mm256_add_epi32(Sum, SumS);

  return _mm256_cvtsi256_si32(Sum);
}

// AVX Quad Accumulators with Pointer Arithmetic
u32 __attribute__((target("avx2"))) QuadAVXPtr(u32 Count, u32 *Input) {
  __m256i SumA = _mm256_setzero_si256();
  __m256i SumB = _mm256_setzero_si256();
  __m256i SumC = _mm256_setzero_si256();
  __m256i SumD = _mm256_setzero_si256();

  Count /= 32;
  while (Count--) {
    SumA = _mm256_add_epi32(SumA, _mm256_loadu_si256((__m256i *)&Input[0]));
    SumB = _mm256_add_epi32(SumB, _mm256_loadu_si256((__m256i *)&Input[8]));
    SumC = _mm256_add_epi32(SumC, _mm256_loadu_si256((__m256i *)&Input[16]));
    SumD = _mm256_add_epi32(SumD, _mm256_loadu_si256((__m256i *)&Input[24]));

    Input += 32;
  }

  __m256i SumAB = _mm256_add_epi32(SumA, SumB);
  __m256i SumCD = _mm256_add_epi32(SumC, SumD);
  __m256i Sum = _mm256_add_epi32(SumAB, SumCD);

  Sum = _mm256_hadd_epi32(Sum, Sum);
  Sum = _mm256_hadd_epi32(Sum, Sum);
  __m256i SumS = _mm256_permute2x128_si256(Sum, Sum, 1 | (1 << 4));
  Sum = _mm256_add_epi32(Sum, SumS);

  return _mm256_cvtsi256_si32(Sum);
}

// ========================================================================
// Infrastructure
// ========================================================================

// Read CPU Timer (RDTSC)
#if defined(__x86_64__) || defined(_M_X64) || defined(__i386__) ||             \
    defined(_M_IX86)
static inline uint64_t read_cpu_timer(void) {
  unsigned int lo, hi;
  __asm__ volatile("rdtsc" : "=a"(lo), "=d"(hi));
  return ((uint64_t)hi << 32) | lo;
}
#else
static inline uint64_t read_cpu_timer(void) { return (uint64_t)clock(); }
#endif

#define ARRAY_SIZE 4096
#define ITERATIONS 100000

#define TEST_FUNCTION(func_name, count, input)                                 \
  do {                                                                         \
    uint64_t min_cycles = UINT64_MAX;                                          \
    volatile u32 result = 0;                                                   \
    for (int iter = 0; iter < ITERATIONS; iter++) {                            \
      uint64_t start = read_cpu_timer();                                       \
      result = func_name(count, input);                                        \
      uint64_t end = read_cpu_timer();                                         \
      uint64_t elapsed = end - start;                                          \
      if (elapsed < min_cycles)                                                \
        min_cycles = elapsed;                                                  \
    }                                                                          \
    double adds_per_cycle = (double)(count) / (double)min_cycles;              \
    printf("%-18s: %8llu cycles, %.6f adds/cycle\n", #func_name,               \
           (unsigned long long)min_cycles, adds_per_cycle);                    \
    fflush(stdout);                                                            \
  } while (0)

int main(void) {
  // Allocate and initialize array
  // Note: malloc guarantees alignment suitable for standard types.
  // For AVX types, using loadu (unaligned load) is safe.
  u32 *input = (u32 *)malloc(ARRAY_SIZE * sizeof(u32));
  if (!input) {
    fprintf(stderr, "Memory allocation failed\n");
    return 1;
  }

  // Fill array with random numbers
  srand(42);
  for (int i = 0; i < ARRAY_SIZE; i++) {
    input[i] = rand() % 100;
  }

  printf("=== IPC & SIMD Demo ===\n");
  printf("Array size: %d integers\n", ARRAY_SIZE);
  printf("Iterations: %d\n\n", ITERATIONS);

  printf("--- Scalar: Single Dependency Chain ---\n");
  TEST_FUNCTION(SingleScalar, ARRAY_SIZE, input);
  TEST_FUNCTION(Unroll2Scalar, ARRAY_SIZE, input);
  TEST_FUNCTION(Unroll4Scalar, ARRAY_SIZE, input);

  printf("\n--- Scalar: Broken Dependency Chain ---\n");
  TEST_FUNCTION(DualScalar, ARRAY_SIZE, input);
  TEST_FUNCTION(QuadScalar, ARRAY_SIZE, input);
  TEST_FUNCTION(QuadScalarPtr, ARRAY_SIZE, input);

  printf("\n--- SIMD Versions (SSE/AVX) ---\n");
  TEST_FUNCTION(SingleSSE, ARRAY_SIZE, input);
  TEST_FUNCTION(SingleAVX, ARRAY_SIZE, input);

  printf("\n--- SIMD: Unrolled & Broken Dependency ---\n");
  TEST_FUNCTION(DualAVX, ARRAY_SIZE, input);
  TEST_FUNCTION(QuadAVX, ARRAY_SIZE, input);
  TEST_FUNCTION(QuadAVXPtr, ARRAY_SIZE, input);

  free(input);
  return 0;
}

测试结果与分析

方法	周期数 (Cycles)	吞吐量 (adds/cycle)	说明
Scalar (Single)	3000	1.37	基准，串行依赖严重
Scalar (Unroll 4x)	2024	2.02	编译器优化试图打破依赖
Scalar (Quad Ptr)	1082	3.79	手动打破依赖 + 指针优化
SIMD (SSE)	860	4.76	128位宽，一次处理4个
SIMD (AVX)	444	9.23	256位宽，一次处理8个
SIMD (Quad AVX Ptr)	320	12.80	AVX + 打破依赖 + 指针优化

结论

从 SingleScalar 的 1.37 到 QuadAVXPtr 的 12.80：

通过结合 SIMD 和打破指令依赖链，我们实现了近 10 倍的性能提升。