Stanford-CS149-并行计算-Lec02-笔记-多核&&超线程&&SIMD

1 简单介绍多核CPU(续上一节课)

多核CPU就是将单个CPU复制很多份, 每一份有自己的ALU、上下文寄存器、乱序单元等等:
multi-core-cpu

使用多核CPU进行并行计算时,简单的思路就是将不存在依赖关系的任务分配到不同的核上执行。

2 并行方案3——SIMD

SIMD是Single Instruction Multiple Data的缩写,即单指令多数据。
simd

  1. SIMD的工作原理

    • 一条指令(Single Instruction)
    • 多个ALU并行处理不同的数据(Multiple Data)
    • 每个ALU执行相同的操作,但处理不同的数据
  2. 示例中的结构

    • 8个ALU(ALU 0-7)
    • 一个取指/解码单元(控制单元)
    • 当执行SIMD指令时:
      • 同一指令广播给所有ALU
      • 每个ALU处理自己的数据部分
      • 8个数据同时被处理

需要说明的是, 这多个ALU存在于同一个处理器中, 所以它们可以共享同一个取指/解码单元。

这里再对ALU和Fetch/Decode组件的关系进行补充:

  1. 指令获取和解码过程

    • Fetch

      • 一次获取一条SIMD指令
      • 同时需要获取对应的多个数据操作数
      • 通常需要较宽的数据总线支持
    • Decode

      • 解析SIMD指令类型
      • 确定数据访问模式
      • 准备数据分发给各个ALU
  2. 数据流动

    1
    2
    3
    4
    5
    Fetch/Decode ──┬→ ALU 0 ← Data[0]
    ├→ ALU 1 ← Data[1]
    ├→ ALU 2 ← Data[2]
    └→ ALU 3 ← Data[3]
    ...

    需要注意的是, 虽然Fetch/Decode是统一的控制单元, 但是实际数据获取可能需要多个周期才能完成, 这也是为什么SIMD处理器的设计需要考虑内存系统的优化。

3 实现SIMDAVX2案例代码

SIMD这一思想对应的指令集就是我们常常在硬件评测类视频或文章中提到的AVXxxx系列指令集。AVX2是目前使用最广泛的SIMD指令集。AVX512支持的SIMD指令集的宽度更大, 但是目前还没有普及, 而且因为发热过大还在Intel最近几代CPU被移除了, 这里的案例代码也是基于AVX2的, 其支持的SIMD指令宽度为256位。

查看个人PC是否支持AVX2指令集的方法:

1
cat /proc/cpuinfo | grep avx2

3.1 中常用的函数总结

3.1.1 浮点运算函数

函数名 说明 示例
单精度浮点运算 (float)
_mm256_add_ps 8个单精度浮点数相加 a + b
_mm256_sub_ps 8个单精度浮点数相减 a - b
_mm256_mul_ps 8个单精度浮点数相乘 a * b
_mm256_div_ps 8个单精度浮点数相除 a / b
双精度浮点运算 (double)
_mm256_add_pd 4个双精度浮点数相加 a + b
_mm256_sub_pd 4个双精度浮点数相减 a - b
_mm256_mul_pd 4个双精度浮点数相乘 a * b
_mm256_div_pd 4个双精度浮点数相除 a / b
数据加载存储
_mm256_load_ps 从对齐内存加载8个float load from memory
_mm256_load_pd 从对齐内存加载4个double load from memory
_mm256_store_ps 存储8个float到对齐内存 store to memory
_mm256_store_pd 存储4个double到对齐内存 store to memory
数据类型转换
_mm256_cvtps_pd float转double (低4个) float -> double
_mm256_cvtpd_ps double转float (4个) double -> float
比较操作
_mm256_cmp_ps 8个float比较 a < b, a == b
_mm256_cmp_pd 4个double比较 a < b, a == b
其他常用操作
_mm256_sqrt_ps 8个float平方根 sqrt(a)
_mm256_sqrt_pd 4个double平方根 sqrt(a)
_mm256_max_ps 8个float最大值 max(a,b)
_mm256_min_ps 8个float最小值 min(a,b)

注意事项:

  1. 函数名中的 ps 表示 packed single (单精度浮点)
  2. 函数名中的 pd 表示 packed double (双精度浮点)
  3. 256 表示使用 256 位寄存器
  4. 单精度运算一次处理 8 个数 (256/32=8)
  5. 双精度运算一次处理 4 个数 (256/64=4)

3.1.2 整型运算函数

函数名 说明 示例
32位整数运算 (int)
_mm256_add_epi32 8个32位整数相加 a + b
_mm256_sub_epi32 8个32位整数相减 a - b
_mm256_mul_epi32 4个32位整数相乘(低位) a * b
_mm256_mullo_epi32 8个32位整数相乘(截断) a * b
16位整数运算 (short)
_mm256_add_epi16 16个16位整数相加 a + b
_mm256_sub_epi16 16个16位整数相减 a - b
_mm256_mullo_epi16 16个16位整数相乘 a * b
8位整数运算 (char)
_mm256_add_epi8 32个8位整数相加 a + b
_mm256_sub_epi8 32个8位整数相减 a - b
数据加载存储
_mm256_load_si256 从对齐内存加载整数 load from memory
_mm256_store_si256 存储整数到对齐内存 store to memory
比较操作
_mm256_cmpgt_epi32 8个32位整数大于比较 a > b
_mm256_cmpeq_epi32 8个32位整数相等比较 a == b
位运算
_mm256_and_si256 按位与 a & b
_mm256_or_si256 按位或 a | b
_mm256_xor_si256 按位异或 a ^ b
移位操作
_mm256_slli_epi32 32位整数左移 a << n
_mm256_srli_epi32 32位整数逻辑右移 a >> n
_mm256_srai_epi32 32位整数算术右移 a >> n

函数名说明:

  1. epi32 表示 32 位有符号整数
  2. epi16 表示 16 位有符号整数
  3. epi8 表示 8 位有符号整数
  4. epu32/epu16/epu8 表示无符号整数版本

一次处理的数量:

  • 32位整数:8个 (256/32=8)
  • 16位整数:16个 (256/16=16)
  • 8位整���:32个 (256/8=32)

这些整数运算在图像处理、音频处理等需要大量整数运算的场景中特别有用。例如:

1
2
3
4
// 8个32位整数同时相加的例子
__m256i a = _mm256_set1_epi32(10); // 设置8个整数都为10
__m256i b = _mm256_set1_epi32(20); // 设置8个整数都为20
__m256i result = _mm256_add_epi32(a, b); // 结果中的8个数都是30

这些函数都需要包含 <immintrin.h> 头文件才能使用。

3.2 OpenMP的并行化技术

我学习avx2的时候顺带学习了OpenMP的并行化技术, 这里也顺带总结一下。
这里阐释一下OpenMP, SIMD, AVX三者的关系和区别:

3.2.1 三种技术的定位

技术 层次 并行类型 特点
OpenMP 线程级 多线程并行 跨核心并行
SIMD 指令级 数据并行 单核内并行
AVX 指令级 数据并行 SIMD的具体实现

3.2.2 使用说明

安装:

1
sudo apt-get install libomp-dev

OpenMP

  • 线程级并行化技术
  • 使用多个 CPU 核心
  • 主要通过 #pragma 指令实现
  • 适用场景:
    1
    2
    3
    4
    #pragma omp parallel for
    for(int i = 0; i < 1000; i++) {
    result[i] = compute(data[i]);
    }

3.2.3 OpenMP关键指令说明

  1. #pragma omp parallel

    • 创建一个并行区域
    • 会创建多个线程,默认为CPU核心数量
    • 并行区域内的代码会被所有线程执行
    • 适用场景:需要执行多个并行操作时
      1
      2
      3
      4
      5
      6
      #pragma omp parallel
      {
      // 这里的代码会被所有线程执行
      int thread_id = omp_get_thread_num(); // 获取线程ID
      // ... 并行执行的代码 ...
      }
  2. #pragma omp parallel for

    • 组合指令,专门用于并行化for循环
    • 自动将循环迭代分配给多个线程
    • 会自动处理循环的任务分配和同步
    • 适用场景:单个循环的并行化
      1
      2
      3
      4
      5
      #pragma omp parallel for
      for(int i = 0; i < n; i++) {
      // 循环的每次迭代会被自动分配给不同线程
      result[i] = compute(data[i]);
      }
  3. 两者的主要区别

    • parallel: 创建一个通用的并行区域,区域内的代码被所有线程执行
    • parallel for: 专门针对循环的并行化,自动进行任务分配,每个迭代只被一个线程执行
  4. 性能考虑

    • parallel 创建线程池的开销较大,适合包含多个并行操作的场景
    • parallel for 针对循环优化,对单个循环并行化更高效

3.2.4 AVX2OpenMP组合使用

可以同时使用这些技术来获得最大性能:

1
2
3
4
5
6
7
#pragma omp parallel for        // OpenMP多线程并行
for(int i = 0; i < n; i += 8) {
__m256 a = _mm256_load_ps(&data[i]); // AVX向量加载
__m256 b = _mm256_load_ps(&other[i]);
__m256 c = _mm256_add_ps(a, b); // AVX向量运算
_mm256_store_ps(&result[i], c); // AVX向量存储
}

3.3 案例代码

这里对官方PPT的sin函数进行了补充, 并添加了OpenMP的并行化比较:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
#include <chrono>
#include <cstdio>
#include <immintrin.h> // Intel 内联函数头文件,包含 SSE/AVX/AVX2 等 SIMD 指令集的内联函数
#include <omp.h> // OpenMP 头文件,用于支持多线程并行编程, #pragma omp parallel 和 #pragma omp for 指令需要该头文件

void sinx_scalar(int N, int terms, float *x, float *y) {
for (int i = 0; i < N; i++) {
float value = x[i];
float numer = x[i] * x[i] * x[i];
int denom = 6; // 3!
int sign = -1;
for (int j = 1; j <= terms; j++) {
value += sign * numer / denom;
numer *= x[i] * x[i];
denom *= (2 * j + 2) * (2 * j + 3);
sign *= -1;
}
y[i] = value;
}
}

void sinx_avx(int N, int terms, float *x, float *y) {
int i;
for (i = 0; i < N - 7; i += 8) {
// 从内存加载 8 个 float 到 256 位向量寄存器
__m256 x_vec = _mm256_loadu_ps(&x[i]);
// 初始化 value 向量为 x_vec
__m256 value = x_vec;
// 计算 x^3 向量
__m256 x_n = _mm256_mul_ps(_mm256_mul_ps(x_vec, x_vec), x_vec);
// 创建一个所有元素都是 6.0f 的向量
__m256 denom = _mm256_set1_ps(6.0f);
// 创建一个所有元素都是 -1.0f 的向量
__m256 sign = _mm256_set1_ps(-1.0f);

for (int j = 1; j <= terms; j++) {
// 计算当前项:sign * numer / denom
__m256 term = _mm256_div_ps(_mm256_mul_ps(sign, x_n), denom);
// 将当前项加到 value 上
value = _mm256_add_ps(value, term);
// 更新 numer:numer *= x^2
x_n = _mm256_mul_ps(x_n, _mm256_mul_ps(x_vec, x_vec));
// 更新 denom:denom *= (2j+2) * (2j+3)
denom = _mm256_mul_ps(denom, _mm256_set1_ps((2 * j + 2) * (2 * j + 3)));
// 翻转 sign 的符号
sign = _mm256_xor_ps(sign, _mm256_set1_ps(-0.0f));
}

// 将计算结果存回内存
_mm256_storeu_ps(&y[i], value);
}

// 处理剩余的元素
for (; i < N; i++) {
float value = x[i];
float numer = x[i] * x[i] * x[i];
int denom = 6; // 3!
int sign = -1;
for (int j = 1; j <= terms; j++) {
value += sign * numer / denom;
numer *= x[i] * x[i];
denom *= (2 * j + 2) * (2 * j + 3);
sign *= -1;
}
y[i] = value;
}
}

void sinx_avx_parallel(int N, int terms, float *x, float *y) {
#pragma omp parallel // 创建一个并行区域,默认创建与 CPU 核心数相同的线程
{
#pragma omp for // 并行化 for 循环:不同的线程会同时处理不同的 i 值区间
// 例如有 4 个线程时的处理方式:
// 线程 0 处理:i = 0, 32, 64, ...
// 线程 1 处理:i = 8, 40, 72, ...
// 线程 2 处理:i = 16, 48, 80, ...
// 线程 3 处理:i = 24, 56, 88, ...
for (int i = 0; i < N - 7; i += 8) {
// 每个线程内部使用 AVX 并行处理 8 个数据
// 实现了两级并行:
// 1 级:OpenMP 多线程并行(不同线程处理不同的 i 区间)
// 2 级:AVX SIMD 并行(每个线程内同时处理 8 个数)
__m256 x_vec = _mm256_loadu_ps(&x[i]); // 加载 8 个连续的 float 数
__m256 value = x_vec;
__m256 x_n = _mm256_mul_ps(_mm256_mul_ps(x_vec, x_vec), x_vec);
__m256 denom = _mm256_set1_ps(6.0f);
__m256 sign = _mm256_set1_ps(-1.0f);

for (int j = 1; j <= terms; j++) {
__m256 term = _mm256_div_ps(_mm256_mul_ps(sign, x_n), denom);
value = _mm256_add_ps(value, term);
x_n = _mm256_mul_ps(x_n, _mm256_mul_ps(x_vec, x_vec));
denom = _mm256_mul_ps(denom, _mm256_set1_ps((2 * j + 2) * (2 * j + 3)));
sign = _mm256_xor_ps(sign, _mm256_set1_ps(-0.0f));
}

_mm256_storeu_ps(&y[i], value); // 存储 8 个结果
}

// 处理剩余元素(不能被 8 整除的部分)
#pragma omp for // 剩余元素也使用多线程并行处理
for (int i = N - (N % 8); i < N; i++) {
// 例如 N=1027 时,最后 3 个元素 (1024,1025,1026) 会被不同线程并行处理
float value = x[i];
float numer = x[i] * x[i] * x[i];
int denom = 6;
int sign = -1;
for (int j = 1; j <= terms; j++) {
value += sign * numer / denom;
numer *= x[i] * x[i];
denom *= (2 * j + 2) * (2 * j + 3);
sign *= -1;
}
y[i] = value;
}
}
}

int main() {
const int N = 1000000;
const int terms = 10;
float *x = new float[N];
float *y_scalar = new float[N];
float *y_avx = new float[N];

// 初始化输入数组
for (int i = 0; i < N; i++) {
x[i] = static_cast<float>(i) / N;
}

// 测量标量版本的性能
auto start = std::chrono::high_resolution_clock::now();
sinx_scalar(N, terms, x, y_scalar);
auto end = std::chrono::high_resolution_clock::now();
auto duration_scalar =
std::chrono::duration_cast<std::chrono::microseconds>(end - start);

// 测量 AVX 版本的性能
start = std::chrono::high_resolution_clock::now();
sinx_avx(N, terms, x, y_avx);
end = std::chrono::high_resolution_clock::now();
auto duration_avx =
std::chrono::duration_cast<std::chrono::microseconds>(end - start);

// 测量并行 AVX 版本的性能
float *y_avx_parallel = new float[N];
start = std::chrono::high_resolution_clock::now();
sinx_avx_parallel(N, terms, x, y_avx_parallel);
end = std::chrono::high_resolution_clock::now();
auto duration_avx_parallel =
std::chrono::duration_cast<std::chrono::microseconds>(end - start);
// 打印结果
printf("标量版本耗时:%ld 微秒\n", duration_scalar.count());
printf("AVX 版本耗时:%ld 微秒\n", duration_avx.count());
printf("并行 AVX 版本耗时:%ld 微秒\n", duration_avx_parallel.count());
printf("加速比:%.2f\n",
static_cast<float>(duration_scalar.count()) / duration_avx.count());
printf("并行 AVX 加速比:%.2f\n",
static_cast<float>(duration_scalar.count()) /
duration_avx_parallel.count());

// 验证结果
for (int i = 0; i < N; i++) {
if (std::abs(y_scalar[i] - y_avx[i]) > 1e-5) {
printf("结果不匹配在索引 %d: 标量 = %f, AVX = %f\n", i, y_scalar[i],
y_avx[i]);
break;
}
}

// 清理内存
delete[] x;
delete[] y_scalar;
delete[] y_avx;

return 0;
}

编译执行:

1
2
g++ -mavx2 -O3 -march=native avx.cpp -o avx
./sinx
编译选项 说明 具体作用
-mavx2 启用 AVX2 指令集 - 允许使用 AVX2 向量指令
- 支持 256 位向量运算
- 启用 mm256* 函数
-O3 最高级别优化 - 函数内联
- 循环优化
- 向量化
- 分支预测优化
- 更积极的指令调度
-march=native 针对本机优化 - 自动检测 CPU 特性
- 使用 CPU 支持的所有指令集
- 生成最优化的机器码

结果:

1
2
3
4
5
6
$ ./avx
标量版本耗时:14692 微秒
AVX 版本耗时:3203 微秒
并行 AVX 版本耗时:2813 微秒
加速比:4.59
并行 AVX 加速比:5.22

4. SIMD问题: 如何处理分支

当并行计算多个通道时, 可能有些通道的计算会由于不同的if判断条件导致不同的计算处理方式, 那么如何处理这种情况?

simd-branch

处理方式: 加上掩码滤除即可:
simd-branch-mask

这里的处理会在第一个作业中使用到, 如果感兴趣可以看下第一个作业: https://github.com/ToniXWD/CS149-asst1

这里引入了一些新新概念:

  1. 指令流一致性(Instruction Stream Coherence)

    • 定义

      • 同一指令序列应用于多个数据元素的程序特性
      • 类似于”同步执行”的概念
    • 重要性

      • 对SIMD处理至关重要
      • 能够充分利用SIMD硬件资源
      • 提高并行效率
    • 多核心情况

      • 不要求指令流一致性
      • 每个核心可以独立执行不同的指令流
      • 更灵活的并行处理方式
  2. 发散执行(Divergent Execution)

    • 定义

      • 指令流不一致的情况
      • 不同数据需要执行不同的指令序列
    • 影响

      • 降低SIMD效率
      • 可能导致串行化
      • 浪费硬件资源
  3. 实际应用示例

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    // 一致性执行(好)
    for(int i = 0; i < n; i++) {
    result[i] = a[i] + b[i]; // 所有数据执行相同操作
    }

    // 发散执行(差)
    for(int i = 0; i < n; i++) {
    if(a[i] > 0)
    result[i] = a[i] + b[i]; // 条件分支导致发散
    else
    result[i] = a[i] - b[i];
    }

5 并行方案4——超线程技术的提出

实际CPU运行时, 由于缓存不命中, 会导致CPU的Stall, 因为CPU太快, 而需要等待慢速的L3 Cache加载数据。

这看起来是一个无解的问题,因为CPU的速度不存在能够与之匹敌的缓存。为了解决这个问题,Intel提出了超线程(Hyper-Threading)技术:

  1. 核心思路

    • 在一个物理核心中维护多个(通常是2个)线程的执行上下文
    • 每个线程都有自己独立的寄存器组和程序计数器
    • 共享核心的执行单元、缓存和其他硬件资源
  2. 工作原理

    • 当一个线程因缓存未命中而停顿时
    • CPU可以快速切换到另一个线程继续执行
    • 切换开销很小,因为硬件维护了完整的线程状态
  3. 优势

    • 提高CPU利用率
    • 隐藏内存访问延迟
    • 不需要完整的上下文切换(与操作系统级线程切换相比)
  4. 限制

    • 线程间共享物理核心的资源
    • 如果两个线程都需要相同的执行资源,性能提升有限
    • 不是所有工作负载都能从超线程中受益

6 小结: 目前学习的3中并行技术

6.1 超标量执行(Superscalar)

  • 定义:在指令流中利用指令级并行性(ILP)
  • 特点
    • 在同一核心内并行处理来自同一指令流的不同指令
    • 并行性由硬件在执行过程中自动发现
  • 实现
    • 动态指令调度
    • 乱序执行
    • 分支预测

6.2 SIMD(单指令多数据)

  • 定义:多个ALU由同一指令控制(在核心内)
  • 优势
    • 对数据并行工作负载高效
    • 在多个ALU上分摊控制成本
  • 实现方式
    • 显式SIMD:由编译器完成向量化
    • 隐式SIMD:由硬件在运行时完成向量化

6.3 多核心(Multi-core)

  • 定义:使用多个处理核心
  • 提供:线程级并行性(Thread-level Parallelism)
  • 特点
    • 每个核心可以同时执行完全不同的指令流
    • 通过软件创建线程来向硬件暴露并行性
  • 实现
    • 通过线程API创建和管理线程
    • 每个核心独立执行不同的线程

6.4 超线程技术

  • 定义:在加载数据时,单核执行另一个不需要IO的控制流
  • 优势
    • 提高ALU利用率
  • 实现方式
    • 在CPU单核内创建多个上下文寄存器

6.5 补充: GPU的结构

gpu-structure

GPU的主要结构可以与CPU进行如下类比:

  1. SM (Streaming Multiprocessor)

    • 类比CPU的一个核心
    • 是GPU的主要计算单元
    • 4090有80个SM单元
  2. CUDA Core

    • 类比CPU中的SIMD ALU
    • 每个SM有128个CUDA Core
    • 专门用于并行执行相同指令的不同数据
    • 与CPU的SIMD ALU类似,都是用于数据并行处理
  3. 层级关系

    1
    2
    3
    4
    5
    6
    7
    8
    9
    GPU ─┬→ SM 0 (类比CPU核心) ─┬→ CUDA Core 0 (类比SIMD ALU)
    │ ├→ CUDA Core 1
    │ └→ ... (共128个)

    ├→ SM 1 ─┬→ CUDA Core 0
    │ ├→ CUDA Core 1
    │ └→ ...

    └→ ... (共80个SM)
  4. 主要区别

    • CPU的SIMD ALU通常一次处理4-8个数据
    • GPU的CUDA Core数量更多,可以同时处理更多数据
    • GPU更适合大规模数据并行计算