算法导论第 3 版之多线程算法(三）

最新推荐文章于 2021-08-16 20:39:13 发布

hoping

最新推荐文章于 2021-08-16 20:39:13 发布

阅读量6.5k

点赞数

文章标签：算法多线程 parallel c 优化 n2

<script type="text/javascript"> </script> 用于学习和交流，欢迎指正。

多线程算法 ( 三）

—— 算法导论第 3 版新增第 27 章

Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, Clifford Stein

邓辉译

原文： http://software.intel.com/sites/products/documentation/cilk/book_chapter.pdf

竞争条件

如果一个多线程算法的计算结果和多核计算机调度指令的方式无关，那么就称其为而确定的 。如果多线程算法的行为在多次运行中会有所不同，那么就称其为非确定的 。一个原本打算成为确定性的多线程算法往往会因为其中包含有“确定性竞争”而不能如愿。

竞争条件是并发的祸根。 Therac-25 放射治疗仪就是一个著名的竞争条件 bug 的受害者，其夺取了 3 个人的生命并使得多人受伤； 2003 年的北美停电问题，也使得 5 千万人无电可用。这些恶性的 bugs 是非常难以发现的。即使在实验室中连续测试数天都不会出现问题的软件，也无法避免其在现场运行时偶尔会出现崩溃。

当两个逻辑上并行的指令去访问同一个内存位置并且只要有其中之一是写入指令时，就会出现确定性竞争 。下面的过程中包含了一个竞争条件：

RACE-EXAMPLE()

1 x = 0

2 parallel for i = 1 to 2

3 x = x+1

4 print x

在第 1 行中，把 x 初始化为 0 ，之后创建了两个并行 strands ，每个都会把 x 增加 1 （第 3 行）。虽然看起来 RACE-EXAMPLE 总是应该打印出值 2 （其串行化版本确实是这样），不过它确实会打印出值 1 。我们来看看为何会出现这种不正常的行为。

当一个处理器在增加 x 时，虽然该操作是不可分割的，但是其是由一系列指令所组成的：

1、从内存中读取 x ，放入处理器的寄存器。

2、把寄存器中的值加 1.

3、把寄存器中的值写回到内存中的 x 。

图 27.5(a) 中展示了 RACE-EXAMPLE 执行过程的计算 dag ，其中的 strands 都被分解成单独的指令。前面讲过，理想的并行计算机支持顺序一致性，因此可以把多线程算法的并行执行看作是满足 dag 中依赖关系的相互交织的指令序列。图中的 (b) 部分展示了一个导致异常行为的计算执行序列中的值。值 x 存储在内存中， r1 和 r2 为寄存器。在步骤 1 中，有个处理器把 x 设置为 0 。在步骤 2 和 3 中，处理器 1 把 x 从内存中读到其寄存器 r1 中，并加 1 ， r1 中的值为 1 。此时，处理器 2 开始执行指令 4 到 6 。处理器 2 把 x 从内存读入到寄存器 r2 ；并加 1 ， r2 中的值为 1 ；接着把该值存入 x ， x 的值为 1 。现在，处理器 1 在步骤 7 时被恢复，把 r1 中的值 1 存入 x ， x 的值其实没有改变。因此，步骤 8 中打印出的是值 1 而不是像串行化版本中的 2 。

我们可以看到所发生的情况。如果并行计算的执行是处理器 1 在处理器 2 执行前执行完了其所有的指令，那么就会打印出值 2 。相反，如果处理器 2 在处理器 1 执行前执行完了其所有的指令，同样会打印出值 2 。但是，当两个处理器的指令同时执行时，就可能像例子中的那样，会丢掉一次对 x 的更新。

当然，有许多种不会导致问题的执行序列。比如，对于像 <1,2,3,4,5,6,7,8> 或者 <1,4,5,6,2,3,7,8> 这样的执行顺序，就会得到正确的结果。这正是确定性竞争的问题所在。通常，大部分的执行顺序会产生正确的结果（比如左边的指令都先于右边的指令的执行，或者相反）。但是当指令交织时，有些顺序会导致不正确的结果。因此，竞争问题非常难以测试。连续数天的测试没有发现任何 bug ，却在现场出现灾难性的系统崩溃（如果后果严重的话）。

虽然应对竞争有多种方法，包括使用互斥锁或者其他同步方法，不过，出于我们的目的，我们只是保证并行运行的 strands 都是独立的 ：它们之间没有任何确定性竞争存在。因此，在 parallel for 中，所有的迭代都应该是独立的。在 spawn 和相应的 sync 之间，被 spawn 的 child 的代码应该和其 parent 的代码独立，包括其他被 spawn 或者被调用的 children 的代码。请注意，传入被 spawn 的 child 的参数是在实际的 spawn 发生之前在其 parent 中求值的，因此被 spawn 的子例程的参数的求值和 spawn 后对这些参数的访问是串行的。

我们用一个例子来说明是多么容易编写出具有竞争问题的代码，下面是一个多线程矩阵 - 向量相乘的有问题的实现，其通过并行化内部 for 循环达成 Θ(lgn) 的 span 值：

MAT-VEC-WRONG(A,x)

1 n = A.rows

2 令 y 为一个新的长度为 n 的向量

3 parallel for i = 1 to n

4 y _i = 0

5 parallel for i = 1 to n

6 parallel for j = 1 to n

7 y _i = y _i + a _ij x _j

8 return y

遗憾的是，这个过程是错误的，因为在第 7 行更新 y _i 时存在竞争问题（对于 j 的 n 个值并行去更新）。练习 27.1-6 会让读者提供一个具有 Θ(lgn) 的 span 值的正确实现。

具有竞争问题的多线程算法有时也是正确的。例如，两个并行线程可以把相同的值存入一个共享的变量中，谁先谁后无关紧要。不过我们通常认为具有竞争问题的代码都是非法的。

来自国际象棋程序的教训

我们以一个真实的故事来结束本小节，该故事发生在世界级多线程国际象棋对弈程序 Socrates[81] 的开发期间（时间做了简略）。该程序的原型是在一个具有 32 个处理器的计算机上进行的，但是最终运行在一个具有 512 个处理器的超级计算机上。某天，开发人员对程序进行了一项优化，针对一项重要的在 32 个处理器上基准测试（ benchmark ），把其运行时间从 T₃₂ = 65 秒减少到 T’₃₂ =40 秒。然而，开发人员使用关于 work 和 span 的性能度量得出结论：在 32 个处理器上更快的优化版本，在 512 个处理器上运行时，实际上比原始版本慢一些。结果，他们放弃了“优化”。

他们的分析方法是这样的。程序的原始版本的 work T₁ = 2048 秒， span T_∞ = 1 秒。如果把不等式 (27.4) 看做是等式， T_P = T₁ /P+T_∞ ，并把它当作在 P 个处理器上的近似运行时间，确实可以得出， T₃₂ =2048/32+1 ＝ 65 。如果做了优化， work 变成 T’₁ =1024 秒， span 变成 T’_∞ = 8 秒。采用同样的近似方法，有 T’₃₂ =1024/32 ＋ 8=40 。

但是，当在 512 个处理器上进行计算运行时间时，这两个版本的相对速度会发生转换。此时， T₅₁₂ =2048/512+1 ＝ 5 秒， T’₅₁₂ =1024/512+8 ＝ 10 秒。在 32 个处理器上能够提速程序的优化版本在 512 个处理器上会使得程序慢 2 倍！优化版本的 span 为 8 ，对于 32 个处理器来说，其不是运行时间中的支配项，但是在 512 个处理器时，就变成了支配项，把更多核的优势化为乌有。

这个故事的寓意在于， work 和 span 是一种好的推断性能的方法，而不是一种好的推断运行时间的方法。

练习（略）

27.2 矩阵相乘的多线程算法

本节将研究矩阵相乘的多线程算法，我们在 4.2 节中研究过这个问题的串行运行时间。我们会介绍基于标准的三重嵌套循环以及基于分治法的矩阵相乘多线程算法。

矩阵相乘多线程算法

我们要研究的第一个算法是直接把 SQUARE-MATRIX-MULTIPLY 过程（第 75 页）中的循环并行化：

P-SQUARE-MATRIX-MULTIPLY(A,B)

1 n = A.rows

2 令 C 为一个新的 n×n 矩阵

3 parallel for i = 1 to n

4 parallel for j = 1 to n

5 c_ij = 0

6 for k = 1 to n

7 c_ij = c_ij + a_ik* b_kj

8 return C

现在来分析这个算法，由于该算法的串行化版本就是 SQUARE-MATRIX-MULTIPLY ，因此其 work 为 T₁ (n)= Θ(n³ ) ，和 S QUARE-MATRIX-MULTIPLY 的运行时间一致。其 span 为 T _∞ (n)= Θ(n) ，因为执行过程先沿着第 3 行启动的 parallel fo r 循环所产生的递归树向下，然后再沿着第 4 行启动的 parallel fo r 循环所属产生的递归树向下，接着执行了第 6 行中普通 for 循环中的所有 n 个迭代，从而总的 span 为， Θ(lgn)+ Θ(lgn) ＋ Θ(n) ＝ Θ(n) 。因此，其 parallelism 为 Θ(n³ )/ Θ(n)= Θ(n² ) 。练习 27.2-3 会要求读者并行化内层循环得到一个 span 为 Θ(lgn) 的算法，不能直接使用 parallel for ，因为会产生竞争问题。

矩阵相乘的多线程分治算法

如 4.2 节中所讲，可以采用 Strassen 分治策略在 Θ(n^lg7 )= Θ(n^2.81 ) 的时间内串行地完成 n×n 矩阵的相乘，这激发我们去寻找该算法的多线程版本。和 4.2 节中一样，我们先来多线程化一个简单一些的分治算法。

回忆一下第 77 页中的 SQUARE-MATRIX-MULTIPLY-RECURSIVE 过程，它把两个 n×n 矩阵 A 、 B 相乘得到一个 n×n 矩阵 C ，采用的方法是把这三个矩阵都分割成四个 n/2×n/2 的子矩阵：

A={{A₁₁ ,A₁₂ },{ A₂₁ ,A₂₂ }}, B={{B₁₁ ,B₁₂ },{ B₂₁ ,B₂₂ }}, C={{C₁₁ ,C₁₂ },{ C₂₁ ,C₂₂ }} 。

我们可以把矩阵积记为：

{{C₁₁ ,C₁₂ },{ C₂₁ ,C₂₂ }} = {{A₁₁ ,A₁₂ },{ A₂₁ ,A₂₂ }} * {{B₁₁ ,B₁₂ },{ B₂₁ ,B₂₂ }}

= {{A₁₁ B₁₁ , A₁₁ B₁₂ },{ A₂₁ B₁₁ , A₂₁ B₁₂ }} +

{{A₁₂ B₂₁ , A₁₂ B₂₂ },{ A₂₂ B₂₁ , A₂₂ B₂₂ }} (27.6)

因此，把 n×n 矩阵的乘法变成 8 个 n/2×n/2 矩阵的乘法操作和一个 n×n 矩阵的加法操作。下面的伪码使用嵌套并行实现了这个分治策略。和 SQUARE-MATRIX-MULTIPLY-RECURSIVE 不同， P-MATRIX-MULTIPLY-RECURSIVE 把输出矩阵作为参数，以避免无关的矩阵创建工作。

P-MATRIX-MULTIPLY-RECURSIVE(C,A,B)

1 n = A.rows

2 if n == 1

3 c ₁₁ = a₁₁ b₁₁

4 else 令 T 为新的 n×n 矩阵

5 把 A 、 B 、 C 和 T 分割成 n/2×n/2 的子矩阵

A₁₁ 、 A₁₂ 、 A₂₁ 、 A₂₂ ； B₁₁ 、 B₁₂ 、 B₂₁ 、 B₂₂ ； C₁₁ 、 C₁₂ 、 C₂₁ 、 C₂₂ ；

以及 T₁₁ 、 T₁₂ 、 T₂₁ 、 T₂₂

6 spawn P -MATRIX-MULTIPLY-RECURSIVE( C_11, A_11, B₁₁ )

7 spawn P -MATRIX-MULTIPLY-RECURSIVE( C_12, A_11, B₁₂ )

8 spawn P -MATRIX-MULTIPLY-RECURSIVE( C_21, A_21, B₁₁ )

9 spawn P -MATRIX-MULTIPLY-RECURSIVE( C_22, A_21, B₁₂ )

10 spawn P -MATRIX-MULTIPLY-RECURSIVE( T_11, A_12, B₂₁ )

11 spawn P -MATRIX-MULTIPLY-RECURSIVE( T_12, A_12, B₂₂ )

12 spawn P -MATRIX-MULTIPLY-RECURSIVE( T_21, A_22, B₂₁ )

13 P -MATRIX-MULTIPLY-RECURSIVE( T_22, A_22, B₂₂ )

14 sync

15 parallel for i = 1 to n

16 parallel for j = 1 to n

17 c_ij = c_ij + t_ij

第 3 行对基本情形进行了处理，也就是进行 1×1 矩阵的相乘。第 4 到 17 行处理了递归情形。在第 4 行中创建了一个临时矩阵 T ，在第 5 行中吧矩阵 A 、 B 、 C 和 T 分割成 n/2×n/2 的子矩阵。（和第 77 页的 SQUARE-MATRIX-MULTIPLY-RECURSIVE 一样，我们忽略如何用索引计算来表示矩阵的子矩阵部分这个小问题）。第 6 行的递归调用把子矩阵 C₁₁ 设置为子矩阵 A₁₁ 和 B₁₁ 乘积，这样 C₁₁ 就等于等式（ 27.6 ）中所形成其和的两项中的第一个。同样，第 7 到 9 行把 C₁₂ 、 C₂₁ 和 C₂₂ 设置为等于等式（ 27.6 ）中形成其和的两项中的第一项。第 10 行把子矩阵 T₁₁ 设置为子矩阵 A₁₂ 和 B₂₁ 乘积，这样 T₁₁ 就等于形成 C₁₁ 和的两项中的第二个。第 11 到 13 行分别把 T₁₂ 、 T₂₁ 和 T₂₂ 设置为形成 C₁₂ 、 C₂₁ 、 C₂₂ 和的两项中的第二个。前 7 个跌鬼调用时 spawn 出来的，最后一个运行在主 strand 中。第 14 行中 sync 语句保证了第 6 到 13 行中的矩阵积的计算都已经完成，然后在第 15 到 17 行中，使用双重嵌套 parallel for 循环把积 T 和 C 相加。

首先，仿照其原始 SQUARE-MATRIX-MULTIPLY-RECURSIVE 过程的串行化运行时间分析方法，来分析一下 P- MATRIX-MULTIPLY-RECURSIVE 过程的 work M₁ (n) 。在递归情形中，分割用时为 Θ(1) ，然后执行 8 个 n/2×n/2 矩阵的递归乘法，最后耗时 Θ(n² ) 进行两个 n×n 矩阵相加。因此，其 work M₁ (n) 的递归等式为：

M₁ (n) = 8 M₁ (n/2) + Θ(n² )

=Θ(n³ ) （根据 master 定理中的情形 1 ）

也就是说，多线程算法的 work 和 4.2 节中的 SQUARE-MATRIX-MULTIPLY （三重嵌套循环）的运行时间完全渐进相同。

现在来确定 P- MATRIX-MULTIPLY-RECURSIVE 的 span M _∞ (n) ，我们首先观察到用于分割的 span 为 Θ(1) ，其被第 15 到 17 行中的双重 parallel for 循环的 span Θ(lgn) 支配。因为 8 个并行递归调用所操作的矩阵大小完全相同，因此任何一个的 span 都是 8 个中最大的那个 span 。所以， P- MATRIX-MULTIPLY-RECURSIVE 的 span M _∞ (n) 的递归等式为：

M_∞ (n) = M_∞ (n/2) + Θ(lgn) (27.7)

这个不符合 master 定理中的任何情形，但是却满足练习 4.6-2 中的条件。根据练习 4.6-2 ，递归式（ 27.2 ）的解为 M_∞ (n) = Θ(lg² n) 。

既然知道了 P- MATRIX-MULTIPLY-RECURSIVE 的 work 和 span ，就可以计算其 parallelism ，为 M₁ (n)/ M_∞ (n) = Θ(n³ /lg² n) ，非常之高。

Strassen 方法的多线程化

我们根据和第 79 页中相同的方法来多线程化 Strassen 算法，仅使用嵌套并行：

1、把输入矩阵 A 、 B 以及输出矩阵按照等式（ 27.6 ）分割成 n/2×n/2 的子矩阵。使用索引计算，该步骤的 wiork 和 span 均为 Θ(1) 。

2、创建 10 个矩阵 S₁ 、 S₂ 、 … 、 S₁₀ ，每个都是 n/2×n/2 的，值为第一步中创建矩阵的和或者差。使用双重嵌套 parallel for 循环，创建者 10 个矩阵的 work 和 span 分别为： Θ(n² ) 和 Θ(lgn) 。

3、使用第 1 步中创建的子矩阵和第 2 步中创建的 10 个矩阵，递归地 spawn 出 7 个计算来计算 7 个 n/2×n/2 的矩阵积 P₁ 、 P₂ 、 … 、 P₇ 。

4、通过对 P_i 矩阵不同组合进行加、减，同样采用双重嵌套 parallel loop 循环计算出所期望的结果矩阵 C 的子矩阵 C₁₁ 、 C₁₂ 、 C₂₁ 、 C₂₂ 。计算着四个子矩阵的 work 和 span 分别为： Θ(n² ) 和 Θ(lgn) 。

现在来分析这个算法，由于其串行化版本和原始的串行算法一样，因此 work 就是串行化的运行时间，也就是 Θ(n^lg7 ) 。对于 P- MATRIX-MULTIPLY-RECURSIVE 来说，我们可以设计出关于 span 的递归等式。在本例中，虽然 7 个递归调用并行执行，但是由于它们操作的矩阵大小完全一样，从而可以像 P- MATRIX-MULTIPLY-RECURSIVE 中一样得到相同的递归式 (27.7) ，其解为 Θ(lg² n) 。因此，多线程 Strassen 算法的 parallelsism 为 Θ(n^lg7 )/ Θ(lg² n) ，虽然比 P- MATRIX-MULTIPLY-RECURSIVE 的 parallelism 稍微低一些,但也很高了。