软件性能工程-多核编程

1. 并行方案

1.1 Pthread

具体参考 CSAPP 的笔记：https://www.caiwen.work/post/csapp-9

创建的线程将由操作系统内核进行调度。使用线程的好处就是可以灵活地进行控制，缺点就是创建线程的开销太大（可以使用线程池解决，即提前创建好若干线程放到线程池中等待使用。我们不断地去复用线程池中的线程，避免了来回创建又销毁的开销），而且需要人为控制，容易出问题。

1.2 Threading Building Blocks

这是由 Intel 主导的一个多任务并行的库。

我个人感觉这个有点像 Rust 的 Tokio。大概思想可以理解为就是他内部有个线程池（线程池中线程的个数一般和 CPU 核心数量差不多，即一个核心上跑一个线程）。

我们把需要并行的任务当作一个“轻量级的协程”。TBB 将在运行时决定把一个任务放到哪个 CPU 的核上去运行。和单纯的线程相比，创建任务的开销很小，并且任务之间上下文的切换的开销也比较小。同时，TBB 使用了 work stealing 的任务调度算法，使得其可以充分利用 CPU 各个核。

大概的语法如下：

值得一提的是，set_ref_count 设置当前任务数量，应设置为要 spawn 出去的子任务数量再加上 1（即自己本身）。然后后面的 spawn_and_wait_for_all 就是根据前面设置的 set_ref_count 来判断是不是已经 wait 掉所有任务了。

1.3 OpenMP

OpenMP 相当于是一个语法层面的扩展，目前很多编译器已经支持。在编译时添加 -fopenmp 即可启用对 OpenMP 的语法解析，并且编译器会链接 OpenMP 的库。

感觉 OpenMP 和 TBB 一样，就是把线程进行了封装，在运行时进行调度。

不过 OpenMP 的一个优点时，适合快速将已有的代码变成并行的，更简洁一点。

2. 并行度

我们可以根据程序的过程得到一个 DAG：

其中一个点连出去多个点，意味着这些点都可以并行地计算。上图标号为 4 的框中，紫色的点连出去另一个紫色的点，表示 spawn 出去一个任务，连出去绿色的点表示 spawn 完任务之后自己继续执行的代码。最后经过 cilk_sync 同步之后会来到蓝色的那个点。

DAG 上的每一个点都被称为一个 strand，表示不存在 spawn、sync、call、return 这些的指令序列。同时，入度为 0 的点称为 initial strand，出度为 0 的点称为 final strand。同时各个边也根据行为有不同的名称：

我们再看一个 DAG：

假定 DAG 上每个点都耗时一个单位时间。

定义 work 表示总工作量，记为 $T_1$，为所有点的数量。也可以理解为是只使用单核计算的耗时。在这个 DAG 中 $T_1=18$。

定义 span 表示 DAG 上的最长点数路径的长度，记为 $T_{\infin}$。也可以理解为假如我们有无穷多个 CPU 核心进行计算的耗时。我们也称 span 为 critical-path length 或是 computational depth。在这个 DAG 中 $T_{\infin}=9$。

令 $P$ 表示 CPU 核心数量，$T_P$ 表示利用 P 个 CPU 核心并行计算所耗费时间，我们有：

Work Law：$T_P\ge \frac{T_1}{P}$。可以认为，我们把全部的 work 都分到所有的核上，这样最少也是 $\frac{T_1}{P}$ 的时间。而考虑到分配的开销，以及可能调度并没有使所有核都一直处于工作状态，所以时间会有所增加。

Span Law：$T_P \ge T_{\infin}$。

定义 speed up 为 $\frac{T_1}{T_P}$。

• 如果 $\frac{T_1}{T_P} < P$，我们就得到了一个 sublinear speedup。
• 如果 $\frac{T_1}{T_P} = P$，我们就得到了一个 linear speedup。
• 如果 $\frac{T_1}{T_P} > P$，我们就得到了一个 superlinear speedup。在我们目前的简单模型中是不会出现的，因为这和 Work Law 是矛盾的。但是实际情况中会出现，这是因为，我们目前隐含的前提是，并行运行和串行运行时，一个任务的耗时是不变的。而实际上可能会出现这样一种情况，比如我们要处理一个很大的数据，直接串行的话可能会使得缓存命中率不高。而如果把数据切分成小块，在不同的 CPU 核心上并行运行时，可能处理的数据大小恰好能放入缓存中，使得缓存命中率增高，因此耗时变短。

定义 parallelism 为我们能获得的最大的 speedup，也称为并行度。根据 Span Law，parallelism 就等于 $\frac{T_1}{T_{\infin}}$。

3. 并行问题

3.1 缓存一致性

CPU 中每个核有自己的 L1 和 L2 缓存，然后又有共享的 L3 缓存。这样一来，就可能出现缓存一致性的问题，比如一个核将某个地址的数据缓存到自己的 L2 缓存中，另一个核又把这个地址的数据修改了，而前者的 L2 缓存此时并没有更新，这就出现了不一致的问题。

为了使得多核之间的缓存保持一致，我们有 MSI 协议，每个缓存行都有如下的三种状态：

• M：当前缓存行中的缓存块被修改了。MSI 协议要求此时不存在其他的核，使得该缓存行处于 M 或是 S 状态。
• S：其他的核中可能也存在这个缓存行，这个缓存行是被共享的。
• I：当前缓存行是无效的。

那么每当一个核要对数据进行修改时，除了修改自己核的缓存中的数据之外，还把当前缓存行状态置为 M，并将其他核中的该缓存行（如果有）状态置为 I。

缓存一致性又带来了新的缓存不命中情形：Sharing miss，而 Sharing miss 又可以分为如下两种：

• True-sharing miss：一个核修改了某个内存之后，把其他核在该内存处的缓存全部失效。此时其他核再访问此处内存就会产生缓存不命中，就是 true-sharing miss。
• False-sharing miss：有可能一个核修改的内存位置和另一个核后续访问的内存位置不同，但是两个位置是一个 cache line。这就使得前者修改后，后者的缓存也失效了。

3.2 Race

determinacy race：假设有两个并行的程序都在访问同一个内存位置，并且至少有一个程序进行了写操作。

data race：假设有两个并行的程序都_无锁地_访问同一个内存位置，并且至少有一个程序进行了写操作。

有时候没有 data race 但会有 determinacy race，如：

benign race：看似没什么影响的 race，如：

cpp
1
2
3
4
struct {
    char a;
    char b;
} x;

如果一个线程修改 x.a，一个线程修改 x.b，这看似是没什么影响的。但是在一些架构或是编译优化行为上，可能会出现将这个结构体看成是 16 位的整体进行更新。这就导致看似只是修改了 x.a ，但其实又修改了 x.b。不过这在 Intel 的 X86 架构上应该是不会出现这种问题的。但是需要注意一点，我们在之前优化技巧中提到了位域语法，使用这个的话就需要注意 Race 的问题了。

3.2.1 Mutex

mutex 有如下属性：

• yielding/spinning
  • yielding：在等待锁的时候会主动向操作系统交出当前进程的控制权
  • spinning：在等待锁的时候会死循环，直到把当前时间片消耗完
• reentrant/nonreentrant
  • reentrant：一个 reentrant mutex 可以做到让一个已经获得到当前锁的线程再次尝试获取当前锁
  • nonreentrant：一个 nonreentrant mutex 会使得已经获得到当前锁的线程再次尝试获取当前锁时会出现死锁
• fair/unfair
  • fair：一个 fair mutex 会使得等待锁实现最长的线程拿到锁
  • unfair：unfair mutex 则不保证这一点

一个 mutex 的实现：

其中 xchg mutex, %eax 是关键。xchg 是一个原子指令，不会被打断，其作用是将 mutex 和 %eax 进行交换。

其实整个 mutex 的关键在 Get_Mutex 中，理论上讲我们并不需要 Spin_Mutex 也能够实现，但是这么做的话我们会不断地执行 xchg 指令修改 mutex，而 mutex 是一个多线程共享的变量，根据 MSI 协议，这会使得其他核上的存有 mutex 的 cache line 被反复失效，引发效率问题。

锁竞争

拿锁时我们希望做到如下的两点：

• 当锁被释放时尽可能拿到锁
• 当锁没被释放时尽可能让浪费的时间少

如果是 spining mutex，会使得锁被释放后立刻拿到锁，但是会使得非常多的时间浪费在拿锁上。如果是 yielding mutex，那么会使得浪费的时间很少，但是当前线程被系统调度出去之后，可能要等很久之后才会再次被调度，锁被释放后不能立即拿到锁。

为了平衡这两点，我们可以有这样一个解决方案：先 spin 一段时间，然后再 yield。spin 的时间与线程被操作系统调度出去再调度回来的时间一致。（lecture 上说是比较优的，我个人感觉和 yielding 差不多）

3.2.2 Reducer

锁太昂贵。在一些情况下我们可以不用加锁，改用一个叫做 Reducer 的东西。

Reducer 的大概思想就是，如果多个线程需要操作同一个全局变量，那么我们把这个全局变量”局部化“，让线程的操作仅自己可见。然后再最后同步的时候再把各个线程对这个全局变量的操作进行合并。但是 Reducer 有两点局限性：首先它要求操作必须有结合律，其次要求线程不需要获取全局变量当前的值。

3.3 Scheduling

现在比较好的一个任务调度策略是 Work Stealing。大概的思想是，每个核心都维护一个 deque，表示当前的任务队列。如果当前任务产生了新的任务，那么就把新任务放入队列的头部，执行这个新任务。如果当前队列的任务都执行完了，那么就去考虑从其他队列的尾部 steal 一些任务过来，使得自己永远保持忙碌。

这个策略的耗时大概是 $T_P \approx \frac{T_1}{P} + O(T_{\infin})$。简单的证明就是，所有核的耗时加起来一定是 $T_1$ 再加上 steal 产生的开销。每次 steal 都有 $\frac{1}{P}$ 的概率 steal 到 critical path 上的任务，使得 span 减一。那么全部的 steal 开销期望就是 $O(PT_{\infin})$。这些时间能分摊到 $P$ 个核心上，所以耗时就是：

$$\frac{T_1+O(PT_{\infin})}{P} = \frac{T_1}{P} + O(T_{\infin})$$

我们也能感性地去考虑为什么这个策略很优。首先，一个基本的思想是，你让所有的核心都在不停地执行任务，这样利用率肯定好。沿着这个思路，大概有两种实现，一种是这里的 work stealing，另一种是，如果有新任务产生了，那么我们考虑把这个新任务分给当前压力最小的队列中去。但后者这么做的话，每次产生一个任务都需要一个分配任务的开销，但 work stealing 的策略只会在出现无事可做之后才会去产生分配的开销。而且后者把任务放到其他队列的时候需要加锁，但 work stealing 可以把锁的开销尽可能减小甚至无锁。（参考论文 The Implementation of the Cilk-5 Multithreaded Language ，链接 https://pages.cs.wisc.edu/~swilson/cilk.pdf 第 7 页。不过论文中似乎只讨论了 Thief 和 Worker 的竞争而没说多个 Thief 抢一个 Worker 该怎么搞）

3.4 堆分配

blowup：在并行时所消耗的堆空间与串行时所消耗堆空间的比值。

user footprint：程序所分配出去的最大空间。

allocator footprint：程序申请到的最大空间。

3.4.1 策略一：Global Heap

直接给 global heap 设置一个锁，当一个线程想要申请/释放内存的时候需要加锁。

这样的 blowup 为 1，但是锁的开销很大。

3.4.2 策略二：Local Heaps

每个线程都有自己的堆。

这样的好处是不需要加锁，没有同步的开销。

坏处是可能发生 memory drift，即在 A 线程申请的内存被 B 线程释放时（比如生产-消费的模式），释放的内存会放入 B 线程中。然后设想一个情况，A 线程一直在申请内存，B 线程一直在释放内存，如果是串行情况的话不会有问题，但并行时发生 memory drift 的话会使得线程 A 的栈空间一直不够用，之前申请的空间全部跑到线程 B 了，从而导致 blowup 没有限制。

3.4.3 策略三：Local Ownership

在策略二的基础上，我们对每个对象记录一下他是从哪个线程申请过来的。释放对象之后把内存放入其所属线程的堆中。

好处是这样 blowup 得到了限制，并且还减少了 false sharing，因为对象都会交给其所属线程进行处理。如果交给不同的线程的话可能会 invalidate 掉整个 cache line。

3.4.4 The Hoard Allocator

有一个 global heap 和 $P$ 个 local heaps。local heap 与 global heap 传递内存的时候以 $S$ 大小的 superblock 为单位。

申请内存时：

不同的是，释放内存的时候，释放掉的内存块会放入当前线程的堆中。当然这么做的话似乎又会出现与策略二一样的问题，因此我们还需要维护如下的式子：

$$u_i \ge min(a_i - 2S, \frac{a_i}{2})$$

其中 $a_i$ 为当前堆拥有的内存大小，$u_i$ 为当前堆实际分配出去的内存大小。

这样做的目的在于，我们可以确保 blowup 为 $O(1+\frac{SP}{U})$。

证明：

$u_i \ge min(a_i - 2S, \frac{a_i}{2})$ 中，$u_i\ge a_i-2S$ 表明一个堆最多有 $2S$ 的空间是没有被分配的，那么 $P$ 个堆最多是 $2SP$ 的空间没有被分配，即 $O(SP)$。$u_i \ge \frac{a_i}{2}$ 表面一个堆最多只能有一般的空间没有被分配，$P$ 个堆最多就是 $O(U)$ 的空间没有被分配。于是有：

$$A \le O(U+SP) \\ \frac{A}{U} = O(1+\frac{SP}{U})$$

我们这个方案相比于方案三，不仅释放内存的时候锁的开销比较少，而且 blowup 还更低。

3.5 内存模型

3.5.1 Relaxed Consistency

处理器发射 store 指令的速度快于 memory network 可以处理的速度。所以在处理器和 memory system 之间会有一个 store buffer。

load 指令拿到数据后，随后的指令会依赖于这个数据，因此 load 指令可能会让处理器暂停，直到从内存中拿到数据。为了让 load 指令的延迟被 cover 掉，现代处理器可能会让 load 指令先于 store 指令执行（这就出现了乱序执行，也就是上图的 load bypass）。如果 load 要读的地址刚好是之前 store 写入的地址，但是 store 操作还位于 store buffer 中的话，load 指令可以直接从 store buffer 中拿到之前写的数据。

我们考虑下面两个线程的代码：

这两个代码可能会交叉执行，但是有一点似乎是应该保证的，那就是一定先执行指令 1 然后再执行指令 2，指令 3 和指令 4 同理。我们称这种保证是 sequential consistency。因此所有可能的情况如下：

因此不可能出现最后 %eax 和 %ebx 都为 0 的情况。

而上面我们讲处理器可能会让 load 先于 store 执行，这就不是 sequential consistency，而是 relaxed consistency。由此可能带来的结果：

Intel X86-64 的内存模型有如下几个规则：（同时可以参考：Manuals for Intel® 64 and IA-32 Architectures 中 Volume 3A 的 10.2 的 Memory Ordering，Intel SDM 之 Memory Ordering - 知乎）

• 同类型操作之间不会发生 reorder
• store 不会被移到 load 的前面（这反而会拖慢处理器运行）
• load 会移动到 store 前面，但前提是 load 和 store 的地址不一样（否则再 reorder 就不对了）
• lock 操作（和锁有关的指令，比如一些原子操作指令）不会与 load 或是 store 进行 reorder

上述规则是对于单个 processor 的。对于多个 processor 之间还有如下的规则：

• 两个 processor 的 store 操作的先后顺序，在不同的 processor 中可能是不一样的，比如下图，在 processor 0 看来，mov [_x], 1 最先执行的，mov [_y,1] 是最后执行的，这导致 r2=0。而在 processor 1 看来，mov [_y,1] 是最先执行的，mov [_x], 1 最后执行的，这导致 r4=0。出现这种问题的原因在于这两个 processor 各自的 store 操作只是被放入了自己的 store buffer 中，还没有被应用到 memory system 中。

• 操作的顺序遵循因果律，比如下图，假设有 r1=1 和 r2=1 成立，那么就说明 processor 0 的 mov [_x], 1 先于 processor 1 的 mov r1, [_x] 执行，然后我们知道 mov r1, [_x] 与 mov [_y], 1 是不能 reorder 的，所以 mov [_y], 1 在随后执行。由于 r2=1，那么同理，后面是执行 mov r2, [_y] 再执行 mov r3, [_x]。于是所有指令的执行先后顺序都被我们推理出来了，并且得知 mov [_x], 1 最先执行，mov r3, [_x] 最后执行，所以 r3 必然为 1.

• 除了执行 store 的处理器之外，所有处理器看到的两个 store 操作的顺序都是一致的。不会出现某些处理器先看到 store a 再 store b，另一些处理器先看到 store b 再看到 store a。如下图，如果假设 r1=1 且 r2=0，那么就说明对于 processor 2 而言，mov [_x], 1 先执行，mov [_y], 1 后执行。那么对于 processor 3 而言也应该是这个顺序，所以就不会出现 r3=1 且 r4=0 的情况。

• lock 操作的顺序对于所有处理器来说都是一样的。有点类似于上面那个规则。

值得一提的是，reorder 不是必须的。如果把 load 提前的话，那么拿到的数据需要在寄存器中存放更久，占用寄存器更多。如果寄存器不够分配了，那么就不会产生 reorder。

3.5.2 Memory Fence

这样的 relaxed consistency 的内存模型会给我们带来麻烦。比如我们实现 peterson 算法（仅用 load 和 store 操作实现锁，见 2025春操作系统训练营游记（一二阶段） - Caiwen的博客）这个算法依赖于 sequential consistency。在 relaxed consistency 下会出现这样的情况：

读取 B_wants 的操作被处理器或是编译器提前到了写 A_wants 之前。对于右侧的代码也是同理。

为了防止代码被 reorder，我们需要 memory fence。memory fence 下方的代码不会被 reorder 到上方。memory fence 的开销比较大，甚至比锁的开销还大（锁的内部应该是有用到 memory fence。而反而是锁的开销更小的原因可能是处理器的设计人员对锁进行了更多的优化）。X86 中有 mfence 指令。在 C 语言中可以使用来自 stdatomic.h 的 atomic_thread_fence() 来实现。

这还没完，不仅处理器会 reorder ，编译器也会，有可能编译器会把我们 critical section 中的操作给 reorder 到上方。我们可以使用 asm volatile("":::"memory") 来实现只针对于编译器的 memory fence：

在 C11 标准中引入了 atomic_store 和 atomic_load。这两个东西保证了读写的原子性（在 32 位系统中对 64 位的读写可能是分两次的），并且确保不会被重排序（ 默认情况，可选内存模型）。被 atomic_store 和 atomic_load 操作的变量需要有 _Atomic 修饰。

4. 并行算法

4.1 循环

考虑如下代码：

cpp
1
2
3
for (int i = 0; i < n; i++) {
    a[i] += b[i];
}

如果我们考虑并行化的话，可以写成如下代码：

其对应的 DAG 为：

我们分析一下：

• Work：总的运算量应该是 $O(n)$ 的
• Span：Span 所在路径为从 Initial Strand 到最右边的灰点，为 $O(G+\frac{n}{G})$

如果我们不使用循环展开，那么 $G=1$，Span 就为 $O(n)$，并行度就是 $O(1)$ 的，约等于没有产生优化。所以如果循环中的 Work 比较小，我们不应该类似于上面这样单纯顺序地循环去 spawn。

如果我们使用循环展开，显然取 $G=\sqrt{n}$ 能够使得 Span 最小，最小值为 $O(\sqrt{n})$，那么并行度为 $O(\sqrt{n})$。

对于循环，更好地做法应该是把循环转化成一种递归的形式：

上面的 G 表示剩余的大小小于 G 就不再继续往下分了。

其对应的 DAG 为：

• Work：$n\cdot I + (\frac{n}{G}-1)\cdot S$ （tips：完全二叉树内部节点数量）
• Span：$G\cdot I + \log {\frac{n}{G}}\cdot S$

粗略地分析，会发现并行度还是比较优的。

4.2 矩阵转置

一般我们是并行化外部循环的，如果只并行化内部循环，span 并不会降低太多。

如果给内层循环也并行化的话，会进一步降低 span：

我们获得的并行度就更大了。一般我们只需要让并行度是 CPU 核心数的 10 倍即可。过大的并行度我们其实并不能充分地利用。反而，由于我们 Spawn 了很多内部节点，所以 Work 有所增加，反而拖慢了我们的速度。

4.3 归并排序

传统归并排序会将当前序列平均分成两部分，分别进行排序，然后再合并。现在我们考虑将两个子部分分别排序的操作并行化。这样做的 Work 是 $O(n\log n)$ 的，Span 为 $T_{\infin}(n) = T_{\infin}(\frac{n}{2})+O(n)=O(n)$，于是并行度为 $O(\log n)$ 的。

Span 的瓶颈在于将两个排好序的序列进行合并。我们考虑也将合并操作并行化：我们令 $A$ 序列为长度较大的那个序列，然后取 $A$ 序列中中间的那个元素，我们假设其为 $x$。那么 $A$ 序列就可以被平均分成大于 $x$ 和小于 $x$ 的两部分。同理 $B$ 序列也可以被不平均分成大于 $x$ 和小于 $x$ 的两部分。其中小于 $x$ 两部分之间的合并和大于 $x$ 两部分之间的合并是可以相互独立，并行执行的。

于是合并操作的 Span 为 $T_{\infin}\le T_{\infin}(\frac{3n}{4})+O(\log n) = O(\log ^2 n)$，Work 为 $O(n)$ （至于怎么来的，我看 ppt 没看太懂...）。

此时归并排序的 Work 不变（合并操作的 Work 都没变），Span 为 $T_{\infin}(n) = T_{\infin}(\frac{n}{2})+O(\log ^2 n) = O(\log ^3 n)$。此时获得的并行度为 $O(\frac{n}{\log ^ 2 n})$。

4.4 二维 DP

假设我们有类似于下面的 DP 转移方程：

$$A_{i,j} = f(A_{i,j-1}, A_{i-1,j}, A_{i-1,j-1})$$

现在我们想将其并行化，可以考虑将整个的 $n\times n$ 区域划分为如下四个部分：

其中同种颜色就是可以并行的。

• Work：$T_1(n)=4T_1{\frac{n}{2}}+O(1)=O(n^2)$
• Span：$T_{\infin}(n)=3T_{\infin}(\frac{n}{2})+O(1)=O(n^{\log 3})$
• 并行度：$O(n^{0.41})$

当然我们还可以做更细的划分：

可以得到 $O(n^{0.53})$ 的并行度。但是这种做法的缓存局部性要差一点，实际中未必比上面那种快。

4.5 BFS

我们把 BFS 的每一层都看作是一个 frontier，我们一个一个 frontier 地处理，同一个 frontier 的点并行处理：

直接并行的话可能会出现 data race，因为 frontier 中不同的点可能会指向同一个点，此时这个点的 parent 设为谁就会有问题。上面的 compare-and-swap 就相当于一个无锁的解决该 race 的方法。

但是这也带来了一个 nondeterministic，对于同一个图进行 BFS，最后得到的 parent 数组都可能不一样。为了再去掉这个 nondeterministic，我们规定如果有多个点都指向同一个点，那么编号最小的点作为该点的 parent。

phase 1 中我们让最小的值才能写入 parent（注意这里不去判断 ngh 是否已经被遍历过）。phase 2 中再遍历一遍所有的边，如果当前点已经作为 ngh 的 parent 之后，我们就把 ngh 的 parent 取负，这样就使得 phase 1 是正确的。

对于大多数图，frontier 的大小呈先增大后减小的趋势：

在 BFS 的时候有两种方式，一种是从 frontier 的点出发，去给别的点尝试标 parent（就是上面的方法），我们称这种是 Top-down 的。另一种是遍历所有还没有设置 parent 的点，然后判断这些点是否有一条连向当前 frontier 的边，有的话就说明为这个点找到了他的 parent，设置并 break，我们称这种是 Bottom-up 的。

如果后者在进行 for all neighbors ngh of v 的时候是按点编号从小到大进行遍历的话，那么这种方式甚至不会产生 nondeterministic，也就不需要原子操作。在 BFS 的早期，遍历到的点还比较少时，使用 Top-down 是优的。随着 BFS 的进行 frontier 可能增大，剩余还没遍历到的点变少，Bottom-up 可能是更优的。所以我们可以结合两种方法，比如当 frontier 大小大于 $\frac{n}{20}$ 的时候使用 Bottom-up，反之则使用 Top-down。