算法导论第3版新增第27章：多线程算法（完整版）

多线程算法（完整版）

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

 

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

 

邓辉 译

 

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

 

       本书中的主要算法都是顺序算法 ，适合于运行在每次只能执行一条指令的单处理器计算机上。在本章中，我们要把算法模型转向并行算法 ，它们可以运行在能够同时执行多条指令的多处理器计算机中。我们将着重探索优雅的动态多线程算法模型，该模型既有助于算法的设计和分析，同时也易于进行高效的实现。

 

       并行计算机（就是具有多个处理单元的计算机）已经变得越来越常见，其在价格和性能方面差距甚大。相对比较便宜的有片上多处理器 桌面电脑和笔记本电脑，其中包含着一个多核集成芯片，容纳着多个处理“核”，每个核都是功能齐全的处理器，可以访问一个公共内存。价格和性能都处于中间的是由多个独立计算机（通常都只是些 PC 级的电脑）组成的集群，通过专用的网络连接在一起。价格最高的是超级计算机，它们常常采用定制的架构和网络以提供最高的性能（每秒执行的指令数）。

 

       多处理器计算机已经以各种形态存在数十年了。计算社团早在计算机科学形成的初期就选定采用随机存取的机器模型来进行串行计算，但是对于并行计算来说，却没有一个公认的模型。这主要是因为供应商无法在并行计算机的架构模型上达成一致。比如，有些并行计算机采用共享内存 ，其中每个处理器都可以直接访问内存的任何位置。而有些并行计算机则使用分布式内存 ，每个处理器的内存都是私有的，要想去访问其他处理器的内存，必须得向其他处理器发送显式的消息。不过，随着多核技术的出现，新的笔记本和桌面电脑目前都成为共享内存的并行计算机，趋势似乎倒向了共享内存多处理这边。虽然一切还是得由时间来证明，不过我们在章中仍将采用共享内存的方法。

 

       对于片上多处理器和其他共享内存并行计算机来说，使用静态线程 是一种常见的编程方法，该方法是一种共享内存“虚拟处理器”或者线程的软件抽象。每个线程维持着自己的程序计数器，可以独立地执行代码。操作系统把线程加载到一个处理器上让其运行，并在其他线程需要运行时将其换下。操作系统允许程序员创建和销毁线程，不过这些操作的开销较大。因此，对于大多数应用来说，在计算期间线程是持久存在的，这也是为何称它们为“静态”的原因。

 

       遗憾的是，在共享内存并行计算机上直接使用静态线程编程非常的困难且易于出错。原因之一是，为了使每个线程所承担的负载大致相当，就需要动态地在线程间分配工作，而这是一项极其复杂的任务。除了那些最简单的应用之外，程序员都得使用复杂的通信协议来实现调度器以对工作进行均衡。这种状况导致了并发平台 的出现，并发平台就是一个用来协调、调度、管理并行计算资源的软件平台。有些并发平台被构建成运行时库，有些则提供了具有编译器和运行时支持的全功能的并行语言。

 

动态多线程编程

       动态多线程 是一种重要的并发平台，也是本章中要采用的模型。使用动态多线程平台，程序员可以无需关心通信协议、负载均衡以及其他静态线程编程中的复杂问题，只要明确应用中的并行性即可。该并发平台中有一个调度器，用来自动均衡计算的负载，因此大大地简化了程序员的工作。动态多线程环境所具有的功能目前还在不断的演化之中，不过基本上都包含有两个功能：嵌套并行以及并行循环。嵌套并行就是可以去“ spawn ”一个子例程，并且使得调用者和子例程能够同时执行。并行循环和普通的 for 循环相似，只是循环中的迭代可以并发地执行。

 

       这两个功能是我们将要在本章中研究的动态多线程模型的基础。该模型的一个关键特征为，程序员只需要指明计算中的逻辑并行性，底层并发平台中的线程会自动调度和均衡计算。我们将研究基于这种模型所编写的多线程算法，以及底层并发平台能够高效进行计算调度的原理。

 

       动态多线程模型具有如下几个重要优点：

l         它是串行编程模型的简单扩展。只需在伪码中增加 3 个“并发”关键字： parallel ， spawn 以及 sync ，就可以描述多线程算法。此外，如果从多线程伪码中去掉这些关键字，就可以得到针对相同问题的串行伪代码，我们称之为“串行化”一个多线程算法。

l         它提供了一种理论上清晰的、基于“ work （工作总量）”和“ span （跨度）”这两个概念量化 parallelism （并行度）的方法。

l         许多多线程算法所涉及的嵌套并行可以从分治范型自然得出。此外，正如串行分治算法可以容易地通过递归关系进行分析一样，多线程算法也是如此。

l         该模型符合并行计算实践的演化方向。越来越多的并发平台开始支持动态多线程技术的不同变种，包括 Cilk [51, 118], Cilk++ [72], OpenMP [60], Task Parallel Library [230], and Threading Building Blocks [292] 。

 

       在 27.1 小节中，我们会介绍动态多线程模型，以及有关 work 、 span 以及 parallelism 的度量方法，我们会使用该度量去分析多线程算法。在 27.2 小节中，我们将研究如何使用多线程来进行矩阵相乘，在 27.3 小节中，我们将处理一个更为困难的问题：归并排序的多线程算法

27.1 动态多线程技术基础

       我们将以递归地计算 Fibonacci 数为例来开始对于动态多线程技术的探索历程。先来回忆一下 3.22 小节中给出的 Fibonacci 数的递归定义：

       F₀ = 0,

       F₁ = 1,

       F_i = F_i-1 + F_i-2  for i ≥ 2.

 

       下面是一个简单的用于计算第 n 个 Fibonacci 数的递归串行算法：

       FIB(n)

       1  if n ≤ 1

2      return n

3  else x = FIB(n-1)

4      y = FIB(n-2)

5      return x+y

 

如果要计算很大的 Fibonacci 数，是不能使用该算法的，因为其中有大量的重复计算。图 27.1 展示了在计算 F₆ 时所创建的递归过程实例树。其中，对于对于 FIB(6) 的调用会递归地调用 FIB(5) 和 FIB(4) 。而对 FIB(5) 的调用又会去调用 FIB(4) 。这两个 FIB(4) 实例返回的结果完全相同（ F₄ =3 ）。由于 FIB 并没有去记住这些结果，因此对于 FIB(4) 的第二次调用重复了第一次调用的工作。

      

       我们用 T(n) 表示 FIB(n) 的运行时间。由于 FIB(n) 包含了两个递归调用和其他一些常数时间的工作，因此得到如下递归方程：

 

       T(n) = T(n-1) + T(n-2) + Θ(1)

 

       我们可以采用替换方法得到该方程的解： T(n) = Θ( Fn) 。作为归纳假设，我们假设 T(n ) ≤aFn-b ，其中a>1 ，b>0 且都为常数。通过替换，我们得到：

 

       T(n) ≤ （aF _n-1 -b ）+ (aF _n-2 -b ）+ Θ(1)

         =   a( F _{n-1 +} F _n-2 ）- 2b + Θ(1)

          =  aF _n -b – (b- Θ(1))

        ≤  aF _n -b

 

    如果我们在选择b 时让其大到足以支配 Θ(1) 中的常量。那么我们接着可以把a 选得大到足以满足初始条件。其分析边界为：

 

    T(n) = Θ( Φⁿ ),                                       (27.1)

 

       其中， Φ=(1+sqrt(5))/2 是黄金分割率，由等式(3.25) 得出。由于F_n 随n 成指数级增长，因此在计算 Fibonacci 数时，该过程非常低效。（问题 31-3 中给出了快得多的方法）。

 

       虽然上面的 FIB 过程对计算 Fibonacci 数来说是一种糟糕的方法，但是在说明多线程算法分析中的关键概念方面，它却是一个好例子。在 FIB(n) 中，第 3 行、第 4 行对 FIB(n-1) 和 FIB(n-2) 的两个递归调用彼此之间相互独立：它们可以按照任意顺序被调用，相互之间也不会有任何影响。因此，这两个递归调用是可以并行运行的。

 

       我们在伪代码中增加了并发关键字 spawn 和 sync 来指示并行属性。下面是采用动态多线程技术重写的 FIB 过程：

P-FIB(n)

       1  if n ≤ 1

2      return n

3  else x = spawn P-FIB(n-1)

4      y = P-FIB(n-2)

5      sync

6      return x+y

 

       请注意，如果我们从 P-FIB 中删除掉并发关键字 spawn 和 sync ，剩下的代码和 FIB 完全一样（除了开始和两处递归调用处的过程名字被更改之外）。我们把多线程算法的串行化 定义为：删除了多线程关键字 spawn 、 sync 以及 parallel （并行循环中会用到）后所得到的串行算法。事实上，我们的多线程伪代码具有一个不错的属性——其串行化版本就是解决相同问题的常用串行伪代码。

 

       在过程调用前面加上 spawn 关键字时，就意味着嵌套并行 ，如第 3 行中所示。 spawn 的语义和普通的过程调用不同，执行 spawn 的过程实例（ parent ）可以和被 spawn 出来的子例程（ child ）并行执行，而不像串行执行中那样去等待 child 执行完成。在本例中，当 child 在计算 P-FIB(n-1) 时， parent 可以并行地去计算第 4 行中的 P-FIB(n-2) 。由于 P-FIB 过程是递归的，因此这两个对其自身进行调用的子例程就创建了嵌套的并行性，对其 children 来说同样如此，于是就产生了一个潜在的巨大子计算树，每个子计算都并行执行。

 

       不过，关键字 spawn 并不是一定要求过程必须和其 child 并发执行，只是表示可以 并发执行。并发关键字表达了计算中的逻辑并行性 ，表明了计算中的哪些部分可以并行的运行。哪些子计算实际上是并发运行的是由调度器 在运行时决定的，在计算进行中，调度器把子计算分配给可用的处理器。稍后，我们会讨论调度器的原理。

 

       一个过程，仅当其执行了 sync 语句时（如第 5 行），才能够安全地使用由其 spawn 的 children 例程的返回值。关键字 sync 表示过程必须等待，直到其所 spawn 的 children 全部完成计算，才能够继续 sync 后面的语句。在 P-FIB 过程中，必须要在 return 语句（第 6 行）前增加 sync 语句，从而避免出现在 x 还没有被计算前就进行 x+y 操作的异常情况。除了 sync 语句所提供的显式同步之外，每个过程都会在其返回前隐式地执行一条 sync 语句，这样可以保证在其终止前，其所有的 children 都已经终止。

 

多线程执行模型

       把多线程计算（由一个代表多线程程序的处理器执行的运行时指令集）看做是一个有向无环图 G= （ V, E ）是很有帮助的，我们称其为计算 dag （ directed acyclic graph ） 。图 27.2

中给出了一个示例，其中的计算 dag 来自计算 P-FIB(4) 。从概念层面来讲， V 中的顶点都是指令， E 中边则表示指令间的依赖关系， (u,v ) ∈ E 表示指令 u 必须在指令 v 之前执行。为了方便起见，如果一条指令链中不包含任何并行控制语句（没有 spawn 、 sync 以及被 spawn 例程中 return ——显式的 return 语句或者过程执行完后的隐式 return ），我们就把它们组成一组，形成一个 strand ，每个 strand 都表示一条或者多条指令。涉及并行控制的指令不包括在 strand 中，但是会出现在 dag 结构中。例如，如果一个 strand 有两个后继，那么其中之一必须得被 spawn 出来，如果一个 strand 有多个前驱，就表示前驱因为一条 sync 语句被合并在一起。因此，一般来说， V 形成了 strand 集合，而有向边集合 E 则表示由并行控制产生的 strand 间的依赖关系。如果 G 中有一个从 strand u 到 strand v 的有向路径，那么我们就说这两个 strands 是（逻辑上）串行的。否则就称其为（逻辑上）并行的。

 

       我们可以把一个多线程计算表示为由内嵌于一棵过程实例树中的 strands 组成的有向无环图。比如，图 27.1 中展示了 P-FIB(6) 的过程实例树，其中没有显示 strands 细节。图 27.2 放大了该树中的一个片段，展现了构成每个过程的 strands 。所有连接 strands 的有向边要么运行于一个过程之中，要么沿着过程树中的无向边运行。

 

       我们可以把计算 dag 的边进行分类，以表示出不同 strands 间依赖关系的种类。在图 27.2 中，沿水平方向连接 strand u 和其同一个过程实例中的后继 u’ 的边被称为 continuation edage （继续边） （ u ， u’ ）。当 strand u spawn 了 strand v 时， dag 中就包含了另一个 spawn edge (u,v) ，在图中显示为指向下的边。表示正常过程调用的 call edge 也指向下。 Strand u spawn 了 strand v 和 u 调用 v 的差别在于： spawn 会产生一条从 u 到其同一过程中后继 u’ 的水平方向的 continuation edge ，意味着 u’ 可以和 v 同时执行，而调用不会产生出这样的边。当 strand u 返回到其调用过程，而 x 是该调用过程中紧跟着下一条 sync 语句的 strand 时，计算 dag 中就会包含 return edge (u,x) ，指向上方。计算从一个 initial strand 开始执行（图 27.2 中被标记为 P-FIB(4) 的过程中的黑色顶点），并以一个 final strand 结束（被标记为 P-FIB(4) 的过程中的白色顶点）。

 

       我们将在理想并行计算机 上研究并行算法的执行，该理想并行计算机由一组处理器和一个顺序一致性 的共享内存组成。顺序一致性的意思是，虽然在实际上多个处理器可以同时对共享内存执行众多的存取操作，但是其产生的结果和在每一步中都只有来自一个处理器的一条指令被执行所产生的完全一样。也就是说，内存的行为就像是按照某个全局的线性顺序来执行指令，该全局顺序保证了每个处理器基于独立的顺序来发出自己的指令。对于动态多线程计算来说，计算是被并发平台自动调度到处理器上的，共享内存的工作方式看起来就像是多线程计算的指令相互交织形成了一个线性的顺序来保持计算 bag 中的偏序关系。这个顺序和调度有关，因此在程序的每次运行可能互不相同，但是每次运行时，我都可以假设指令是按照和计算 bag 一致的某个线性顺序执行的，并基于这个假设来理解其行为。

 

       除了对语义进行假设外，还可以对理想并行计算机模型做一些性能方面的假设。特别地，我们假设机器中的每个处理器具有相同的计算能力，并忽略掉调度的开销。虽然后面这个假设听起来过于乐观，不过在实践中，对于具有充分“ parallelism （并行度）”（后面会准确定义这个术语）的算法来说，调度的开销通常是极其小的。

性能度量

       我们可以使用两个度量：“ work ”和“ span ”，来衡量多线程算法的理论效率。 work 指的是在一个处理器上完成全部的计算所需要的总时间。也就是说， work 是所有 strand 执行时间的总和。如果计算 dag 中每个 strand 都花费单位时间，那么其 work 就是 dag 中顶点的数目。 span 是在沿 dag 中任意路径执行 strand 所花费的最长时间。同样，如果 dag 中每个 strand 都花费单位时间，那么其 span 就等于 dag 中最长路径（也就是关键路径 ）上顶点的数目。（在 24.2 节中讲过，可以在 Θ (V+E) 时间内找到 dag G=(V,E) 的一条关键路径 ）。例如，图 27.2 中的计算 dag 共有 17 个顶点，其中 8 个在关键路径上，因此，如果每个 strand 花费单位时间的话，那么其 work 是 17 个单位时间，其 span 为 8 个单位时间。

 

       多线程计算的实际运行时间不仅依赖于其 work 和 span ，还和可用处理器的数目以及调度器向处理器分配 strand 的策略有关。我们用下标 P 来表示一个在 P 个处理器上的多线程计算的运行时间。比如，我们用 T _P 来表示算法在 P 个处理器上的运行时间。 work 就是在一个处理器上的运行时间，也就是 T ₁ 。 span 就是每个 strand 具有自己独立处理器时的运行时间（也就是说，如果可用的处理器数目是无限的），用 T _∞ 来表示。

 

work 和 span 提供了在 P 个处理器上运行的多线程计算花费时间 T _P 的下界：

l         在一个单位时间中，具有 P 个处理器的理想并行计算机最多能够完成 P 个单位工作，因此在 T _P 时间内，能够完成最多 PT _P 数量的工作。由于总的工作为 T ₁ ，因此我们有： PT _P ≥ T ₁ 。两边同除以P 得到work 法则（work law ） ：

 

T _P ≥ T ₁ /P.                                                           (27.2)

 

l         具有 P 个处理器的理想并行计算机肯定无法快过具有无限数量处理器的机器。换种说法，具有无限数量处理器的机器可以通过仅使用 P 个处理器的方法来仿真具有 P 个处理器的机器。因此，得到 span 法则（ spaw law ） ：

 

       T _P ≥ T _∞ .                                                              (27.3)

 

       我们用比率 T ₁ / T _P 来定义在 P 个处理器上一个计算的加速因子（ speedup ） ，它表示该计算在 P 个处理器上比在 1 个处理器上快多少倍。根据 work 法则， T _P ≥ T ₁ /P ，意味着 T ₁ /T _P ≤ P 。因此，在 P 个处理器上的加速因子最多为 P 。当加速因子和处理器的数目成线性关系时，也就是说，当 T ₁ /T _P = Θ P 时，该计算具有线性加速 的性质，当 T ₁ /T _P =P 时，称其为完全的线性加速 。

 

       我们把 work 和 span 的比率 T ₁ /T _∞ 定义为多线程计算的 parallelism （并行度） 。可以从三个角度来理解 parallelism 。作为一个比率， parallelism 表示了对于关键路径上的每一步，能够并行执行的平均工作量。作为一个上限， parallelism 给出了在具有任何数量处理器的机器上，能达到的最大可能加速。最后，也是最重要的，在达成完全线性加速的可能性上， parallelism 提供了一个在限制。具体地说，就是一旦处理器的数目超过了 parallelism ，那么计算就不可能达成完全线性加速。为了说明最后一点，我们假设 P > T ₁ /T _∞ ，根据 span 法则，加速因子满足 T ₁ /T _P ≤ T ₁ /T _∞ <P 。此外，如果理想并行计算机的处理器数目 P 大大超过了 parallelism （也就是说，如果 P >> T ₁ /T _∞ ），那么 T ₁ /T _P <<P ，这样，加速因子就远小于处理器的数目。换句话说，处理器的数目超过 parallelism 越多，就越无法达成完全加速。

 

       例如，我们来看看图 27.2 中 P-FIB(4) 的计算过程，并假设每个 strand 花费单位时间。由于 work T ₁ =17 ， span T _∞ =8 ，因此 parallelism T ₁ /T _∞ =17/8=2.125 。从而，无论我们用多少处理器来执行该计算，都无法获得 2 倍以上的加速因子。不过，对于更大一些的输入来说， P-FIB(n) 会呈现出更大的 parallelism 。

 

       我们把在一台具有 P 个处理器的理想并行计算机上执行多线程算法的并行 slackness （闲置因子） 定义为： (T ₁ /T _∞ )/P = T ₁ /(PT _∞ ) ，也就是计算的 parallelism 超过机器处理器数目的倍数因子。因此，如果 slackness 小于1 ，那么就不能达成完全的线性加速，因为 T ₁ /(PT _∞ )<1 ，根据 span 法则，在 P 个处理器上的加速因子满足 T ₁ /T _P ≤ T ₁ /T _∞ <P 。事实上，随着 slackness 从 1 降低到 0 ，计算的加速因子就越来越远离完全线性加速。如果 slackness 大于 1 ，那么单个处理器上工作量就成为限制约束。我们将看到，随着 slackness 从 1 开始增加，一个好的调度器可以越来越接近于完全线性加速。

调度

       好的性能并不仅仅来自于对 work 和 span 的最小化，还必须能够高效地把 strands 调度到并行计算机的处理器上。我们的多线程编程模型中没有提供指定哪些 strands 运行在哪些处理器上的方法。而是依赖于并发平台的调度器来把动态展开的计算映射到单独的处理器上。事实上，调度器只把 strands 映射到静态线程，由操作系统来把线程调度到处理器上，不过这个额外的间接层次并不是理解调度原理所必需的。我们可以就认为是由并发平台的调度器直接把 strands 映射到处理器的。

 

       多线程调度器必须能够在事先不知道 strands 何时被 spawn 以及何时完成的情况下进行计算的调度——它必须在线（ on-line ） 操作。此外，一个好的调度器是以分散的（ distributed ）形式运转的，其中实现调度器的线程互相协作以均衡计算负载。好的在线、分散式调度器确实存在，不过对它们进行分析是非常困难的。

 

       因此，为了简化分析工作，我们将研究一个在线、集中式（ centralized ） 调度器，在任意时刻，它都知道计算的全局状态。我们将特别分析贪婪式调度器 ，它们会在每个执行步骤中把尽可能多的 strands 分配给处理器。如果在一个执行步骤中有至少 P 个 strands 可以执行，那么就称这个步骤为完全步骤 ，贪婪调度器会把就绪 strands 中的任意 P 个分配给处理器。否则，如果就绪的 strands 少于 P 个，则称这个步骤为不完全步骤 ，调度器会把每个 strand 分配给独立的处理器。

 

       根据 work 法则，在 P 个处理器上可以达到的最快运行时间为 T _P = T ₁ /P ，根据 span 法则，最好的情况是 T _P = T _∞ 。下面的定理表明，因为贪婪式调度器可以以这两个下界之和为其上界，所以其可被证明是一个好的调度器。

 

定理 27.1

       在一台具有 P 个处理器的理想并行计算机上，对于一个 wrok 为 T ₁ ， span 为 T _∞ 的多线程计算，贪婪调度器执行该计算的时间为：

 

T _P ≤ T ₁ / P + T _∞ .                                                              (27.4)

 

证明： 首先来考虑完全步骤。在每个完全步骤中， P 个处理器完成的工作总量为 P 。我们采用反证法，假设完全步骤的数目严格大于 └ T ₁ / P ┘，那么完全步骤所完成的工作总量至少为：

 

P*( └ T ₁ / P ┘+1) = P └ T ₁ / P ┘ + Ｐ

　　　　　　 　 = 　 T ₁ －（ T ₁ mod P ） + P  ( 根据等式 3.8 得出 )

               > T ₁                               ( 根据不等式 3.9 得出 ) 。

 

因此， P 个处理器所完成的工作比所需要的还多，矛盾，所以完全步骤的数目最多为 └ T ₁ / P ┘。

 

    现在，考虑一个不完全步骤。我们用 G 来表示整个计算的dag ，不失一般性，假设每个strand 都花费单位时间。（我们可以把超过单位时间的strand 用一串单位时间strand 来替代）。令G’ 为在该不完全步骤开始时G 已经执行的部分构成的子图，令G” 为在该不完全步骤完成后G 中还没有执行的部分构成的子图。dag 中最长的路径一定起始于入度（ in-degree ）为0 的顶点。由于贪婪调度器中的一个不完全步骤会把G’ 中所有入度为0 的 strands 全部执行，因此 G ”的最长路径长度一定不 G ’中的最长路径小 1 。换句话说，一个不完全步骤会把还没有执行的 dag 的 span 减 1 。所以，非完全步骤的数目最多为 T _∞ 。

 

       由于每个步骤要么是完全的，要么是不完全的，因此定理得证。

 

       下面是定理 27.1 的推论，说明了贪婪式调度器总是具有好的调度性能。

 

推论 27.2

在一台具有 P 个处理器的理想并行计算机上，任何由贪婪式调度器调度的多线程计算的运行时间 T _P ，不会超过最优时间的 2 倍。

 

证明： 令 T _P ^* 为在具有 P 个处理器的机器上，一个最优调度器产生的运行时间，令 T ₁ 和 T _∞ 为该计算的 work 和 span 。根据 work 法则和 span 法则（不等式 27.2 和 27.3 ），得出：

T _P ^* ≥ max(T ₁ /P, T _∞ ) ，根据定理 27.1 ，有：

T _P ≤ T ₁ /P + T _∞

≤ 2*max(T ₁ /P, T _∞ )

≤ 2* T _P ^*

 

              下一个推论告诉我们，对于任何多线程计算来所，随着 slackness 的增长，贪婪式调度器都可以达到接近完全的线性加速。

 

推论 27.3

令 T _P 为在一台具有 P 个处理器的理想并行计算机上，贪婪式调度器调度一个多线程计算的运行时间，令 T ₁ 和 T _∞ 为该计算的 work 和 span 。那么如果 P << T ₁ /T _∞ ，就有 T _P ≈ T ₁ /P （或者相等），也就是具有大约为 P 的加速因子。

 

证明： 假设 P<< T ₁ /T _∞ ，那么就有 T _∞ << T ₁ /P ，因此根据定理 27.1 ，有 T _P ≤ T ₁ /P + T _∞ 。根据 work 法则（ 27.2 ）得到 T _P ≥ T ₁ /P ，因此得出 T _P ≈ T ₁ /P （或者相等），加速因子为： T ₁ /T _P ≈ P 。

 

       符号 << 表示“远小于”，但是“远小于”意味着多少呢？作为经验之谈，当 slackness 至少为 10 时（也就是说， parallelism 是处理器数目的 10 倍），通常就足以得到很高的加速因子。贪婪调度器的上界不等式（ 27.4 ）中的 span 项小于单处理器 work 项的 10% ，这对于绝大多数实际应用情况而言已经足够好了。例如，如果一个计算仅在 10 个或者 1000 个处理器上运行，那么去说 1,000,000 的 parallelism 比 10,000 更好是没有意义的，即使它们之间有 100 倍的差异。正如问题 27-2 所表明的那样，有时通过降低计算的最大并行度，所得到的算法要好于关注其他问题所得到算法，并且还能在相当数目的处理器上伸缩良好。

多线程算法分析

       现在，我们已经拥有了分析多线程算法的所有工具，并且对于在不同数目处理器上的运行时间也有了个不错的边界。对于 work 的分析相对简单，因为只不过就是分析一个普通的串行算法的运行时间（也就是多线程算法的串行化版本），对此，我们早已熟悉，这正是本书大部分内容所讲的东西！对 span 的分析会更有趣一些，一旦掌握了其中的诀窍，通常也不难。我们将以 P-FIB 程序为例来研究一些基本概念。

       分析 P-FIB(n) 的 work T ₁ (n) 没什么难度，因为我们已经做过了。原始的 FIB 过程就是 P-FIB 的串行化版本，因此 T ₁ (n)= T(n)= Θ( Φⁿ ) （基于等式 27.1 ）。

 

    图27.3 中展示了如何去分析span 。如果两个子计算被串行合并在一起，那么其组合的span 等于二者span 之和，如果它们被并行合并在一起，那么其组合的span 等于二者span 中较大的那一个。对于 P-FIB(n) 来说，第3 行中spawn 的P-FIB(n-1) 和第4 行中spawn 的P-FIB(n-2) 并行运行。因此，我们可以把P-FIB(n) 的span 表示为如下递归式：

T _∞ (n) = max(T _∞ (n-1), T _∞ (n-2)) + Θ(1)

     = T _∞ (n-1) + Θ(1),

结果为： T _∞ (n) = Θ(n) 。

 

    P-FIB(n) 的 parallelism 为 T ₁ (n)/ T _∞ (n) = Θ( Φⁿ /n) ，其随着n 增长的速度极快。因此，对 P-FIB(n) 来说， 即使在最大的并行计算机上，一个中等大小的n 值就足以获得接近完全的线性加速，因为该过程具有相当大的并行 slackness 。

并行循环

       有许多算法，其包含的循环中的所有迭代都是可以并行执行的。我们将看到，可以使用 spawn 和 sync 关键字来并行化这种循环，不过如果能够直接指明这种循环的迭代可以并发执行的话，会更加方便一些。我们通过使用 parallel 并发关键字来在伪码中提供该功能，它位于 for 循环语句的 for 关键字之前。

 

       我们以一个 n ×n 的矩阵A= （a _ij ）乘以一个n 元向量x= （x _j ）为例进行说明。相乘的结果为一个n 元向量y= （y _i ），如下：

 

y _i = ∑ ⁿ _j=1 a _ij x _j ，

 

i=1 ，2 ，…,n 。我们可以通过并行地计算y 的所有项来进行矩阵- 向量的乘法操作，如下：

MAT-VEC(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             for j = 1 to n

7                y _i = y _i + a _ij x _j

8         return y

 

       在这段代码中，第 3 行和第 5 行中的 parallel for 关键字表示这两个循环中的迭代都可以并发执行。编译器可以把 parallel for 循环实现为基于嵌套并行的分治式子例程。例如，第 5 到 7 行中的 parallel for 循环可以被实现为对 MAT-VEC-MAIN-LOOP(A,x,y,n,l,n) 的调用，子例程 MAT-VEC-MAIN-LOOP 是编译器生成的辅助子例程，如下：

MAT-VEC-MAIN-LOOP(A, x, y, n, i, i’)

1               if i == i’

2                    for j = 1 to n

3                        y _i = y _i + a _ij x _j

4               else mid = └(i+i’)/2 ┘

5                    spawn MAT-VEC-MAIN-LOOP(A, x, y, n, i, mid)

6                    MAT-VEC-MAIN-LOOP(A, x, y, n, mid+1, i’)

7                    sync

 

       该代码递归地 spawn 循环中的前半部分迭代，使其和后半部分迭代并行执行，然后执行一条 sync 语句，创建了一棵二叉树式的执行过程，其中叶子为单独的循环迭代，如图 27.4 所示。

       现在来计算对于 n ×n 矩阵， MAT-VEC 的 work T ₁ (n) ，也就是计算其串行化版本的运行时间，这个串行化版本可以通过把 parallel for 循环替换成普通的 for 循环得到。由此，我们得到 T ₁ (n)= Θ(n² ) ，因为第5 到7 行的两重嵌套循环所产生的平方级运行时间占支配地位。在这个分析中，我们忽略掉了实现并行循环的递归 spawn 的开销。事实上，和其串行化版本相比，递归spawn 的开销确实增加了并行循环的工作量，不过并不是渐进关系的。原因如下，因为递归过程实例树是一颗满二叉树，所以内部节点的个数正好比叶子的个数少 1 （见练习 B.5-3 ）。每个内部节点分割迭代范围时所耗费的都是常数时间，并且每个叶子都对应循环中的一个迭代，其至少耗费常数时间（在本例中是Θ(n) ）。因此，我们可以把递归 spawn 的开销分摊到迭代的工作中，对全部工作来说，至多增加了一个常数倍数因此。

 

       在实际实现中，动态多线程并发平台时常会在一个叶子中执行多个迭代，从而使得递归产生的叶子的粒度变粗，这个过程可以是自动地，也可以由程序员来控制，因此减少了递归 spawn 的开销。付出的代价是降低了并行度，不过，如果计算具有局够大的并行 slackness ，那么还是可以达成接近完全的线性加速的。

 

       在分析并行循环的 span 时，也必须得考虑到递归 spawn 的开销。由于递归调用的深度和迭代的次数成对数关系，因此对于一个具有 n 次迭代，其第 i 个迭代的 span 为 iter _∞ (i) 的并行循环来说，其 span 为：

T _∞ (n)= Θ(lgn) + max_{1 ≤i ≤n} iter _∞ (i) 。

 

例如，对于以一个 n ×n 矩阵为参数的 MAT-VEC 来说，第 3 、 4 行中的并行初始化循环的 span 为 Θ (lgn) ，因为和每个迭代中的常数工作时间相比，递归spawn 占支配地位。第5 到7 行中的双重嵌套循环的span 为Θ (n) ，因为外层 parallel for 循环的每个迭代都包含着内层（串行） for 循环的 n 个迭代。伪码中剩余部分的 span 为常数，因此整个过程的 span 由双重嵌套循环支配，也就是 Θ (n) 。由于过程的 work 为 Θ(n² ) ，所以 parallelism 为 Θ(n² )/ Θ(n)  =Θ(n) 。（练习 27.1-6 会让读者提供一个具有更高并行度的实现）。

竞争条件

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

 

       竞争条件是并发的祸根。 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 稍微低一些，但 也非常高了。

练习（略）

27.3 多线程归并排序算法

       我们最先是在 2.3.1 节中见到串行归并排序的，在 2.3.2 节中我们分析了它的运行时间，得出的结果是 n Θ(lgn) 。因为归并排序已经使用了分治范型，所以看起来非常适合采用嵌套并行对其进行多线程化。可以很容易地修改伪码，使第一个递归调用被 spawn 出来：

 

MERGE-SORT’(A, p, r)

1                if p < r

2                     q = └ (p+r) ┘ /2

3                     spawn MERGE-SORT’(A, p, q)

4                     MERGE-SORT’(A, q+1, r)

5                     Sync

6                     MERGE(A, p, q, r)

 

和其串行版本了一样， MERGE-SORT’ 对子数组 A[p..r] 进行排序。在第 3 、 4 行中的两个递归子例程完成后（由第 5 行中的 sync 语句保证）， MERGE-SORT’ 调用了同一个 MERGE 过程（第 31 页）。

 

       我们来分析一下 MERGE-SORT’ 。首先，来分析一下 MERGE 。以前讲过，归并 n 个元素的串行运行时间为 Θ(n) 。因为 MERGE 为串行执行的，所以其 work 和 span 均为 Θ(n) 。因此，如下的递归式刻画了 MERGE-SORT’ 归并 n 个元素时的 work MS’₁ (n) ：

 

MS’₁ (n) = 2MS’₁ (n/2) + Θ(n) = 　 Θ(nlgn) ，

 

和归并排序的串行执行时间一样。由于对 MERGE-SORT’ 的两个递归调用可以并行运行，因此 span MS _∞ ’(n) 由如下递归式定义：

 

MS _∞ ’(n) = MS _∞ ’(n/2) + Θ(n) 　＝　 Θ(n) 。

 

因此， MERGE-SORT’ 的 parallelism 为 MS’₁ (n)/ MS _∞ ’(n)= Θ(lgn) ，不是很令人满意。要对 1 千万个元素进行排序，在处理器不多时，还可能能够获取线性加速，但是无法高效地扩展到数百个处理器的情形。

 

       读者也许已经知道多线程归并排序中 parallelism 的瓶颈所在：串行的 MERGE 过程。虽然初看起来归并好像天生就是串行的，但是事实上，我们可以通过使用嵌套并行将其变成多线程的。

 

       图 27.6 中说明了我们的多线程归并的分治策略，其操作的对象是数组 T 的子数组。假设我们要把两个已排序子数组——长度为 n1=r1-p1+1 的数组 T[p1..r1] 和长度为 n2=r2-p2+1 的数组 T[p2..r2] ——归并成另外一个长度为 n3=r3-p3+1=n1+n2 的数组 A[p3..r3] 。不失一般性，我们假设 n1>n2 。

      

       我们先来找到子数组 T[p1..r1] 的中间元素 x=T[q1] ，其中 q1= └ (p1+r1) ┘ /2 。因为子数组是排过序的，所以 x 是 T[p1..r1] 的中间值：所有 T[p1..q1-1] 中的元素都不大于 x ，所有 T[q1+1..r1] 中的元素都不小于 x 。然后，然后我们使用二分查找找到子数组 T[p2..r2] 中的索引 q2 ，使得把 x 插入到 T[q2-1] 和 T[q2] 之间后，该子数组仍然是有序的。

 

       接下来，我们按照如下步骤来把原始子数组 T[p1..r1] 和 T[p2..r2] 归并为 A[p3..r3] ：

1、 令 q3=p3+(q1-p1)+(q2-p2)

2、 把 x 复制到 A[q3]

3、 递归地归并 T[p1..q1-1] 和 T[p2..q2-1] ，并把结果存入子数组 A[p3..q3-1]

4、 递归地归并 T[q1+1..r1] 和 T[q2..r2] ，并把结果存入子数组 A[q3+1..r3]

 

  在计算 q3 时， (q1-p1) 是子数组 T[p1..q1-1] 中元素的数量， (q2-p2) 是子数组 T[p2..q2-1] 中元素的数量。因此，其和就是子数组 A[p3..r3] 中 x 之前的元素的个数。

 

       n1=n2=0 为基本情形处理，对空的子数组来说，无需任何归并工作要做。由于我们假设子数组 T[p1..r1] 至少和 T[p2..r2] 一样长，也就是说， n 1 ≥ n2 ，所以在基本情形中只需检查是否 n1=0 即可。我们必须得保证，在递归中正确地处理了两个子数组中只有一个为空的情况，根据我们 n 1 ≥ n2 的假设，这个子数组必定是 T[p2..r2] 。

 

       现在，我们来把这些想法变成伪代码。先编写二分查找，以串行方式实现。过程 BINARY-SEARCH(x, T, p, r) 接收一个键值 x 和子数组 T[p..r] ，返回值为下面之一：

l         如果 T[p..r] 为空（ r<p ），那么返回索引 p 。

l         如果 x ≤ T[p] ，因此也小于等于 T[p..r] 中的所用元素，返回索引 p 。

l         如果 x> T[p] ，那么就返回区间 p<q ≤ r+1 中满足 T[q-1]<x 的最大索引 q 。

 

伪代码如下：

BINARY-SEARCH(x, T, p, r)

1                low = p

2                high = max(p, r+1)

3                while low < high

4                     mid = └(low+high)/2┘

5                     if x ≤ T[mid]

6                         high = mid

7                     else low = mid + 1

8                return high

 

BINARY-SEARCH(x, T, p, r) 在最坏情况下的时间复杂度为 Θ(lgn) ，其中 n=r-p+1 是所操作的子数组的长度。（请参见练习 2.3-5 ）。由于 BINARY-SEARCH 是一个串行过程，因此其最坏情况下的 work 和 span 均为 Θ(lgn) 。

 

       现在，我们来编写多线程归并过程的伪码。和第 31 页中的 MERGE 过程一样， P-MERGE 过程同样也假设要归并的两个子数组位于同一个数组中。不过，和 MERGE 不同的是， P-MERGE 不需要假定要归并的两个子数组必须是相邻的。（也就是说， P-MERGE 不要求 p2=r1+1 ）。二者之间另外一个区别是， P-MERGE 多了一个输出参数：子数组 A ，用于存放归并后的值。 P-MERGE(T,p1,r1,p2,r2,A,p3) 把有序数组 T[p1..r1] 和 T[p2..r2] 归并为子数组 A[p3..r3] ，其中 r3 并没有作为参数传入，其值为： p3+(r1-p1+1)+(r2-p2+1)-1=p3+(r1-p1)+(r2-p2)+1 。

 

P-MERGE(T,p1,r1,p2,r2,A,p3)

1  n1 = r1-p1+1

2  n2 = r2-p2+1

3  if n1 < n2    // 确保 n1 ≥ n2

4      交换 p1 和 p2

5      交换 r1 和 r2

6       交换 n1 和 n2

7  if n1 == 0   // 均为空？

8      return

9  else q1 = └(p1+r1)/2┘

10     q2 = BINARY-SEARCH(T[q1],T,p2,r2)

11     q3 = p3 + (q1-p1) +(q2-p2)

12     A [q3]=T[q1]

13     spawn P-MERGE(T,p1,q1-1,p2,q2-1,A,p3)

14     P-MERGE(T,q1+1,r1,q2,r2,A,q3+1)

15     sync

 

       P-MERGE 的工作过程如下。第 1 、 2 行分别计算子数组 T[p1..r1] 和 T[p2..r2] 的长度。第 3 到 6 行确保 n 1 ≥ n2 。第 7 行测试基础情况，也就是子数组 T[p1..r1] 为空（此时 T[p2..r2] 也为空），此时直接返回。第 9 到 15 行实现了分治策略。第 9 行计算 T[p1..r1] 的中点，第 10 行查找 T[p2.r2] 中索引 q2 ，使得 T[p2. .q2-1] 中的所有元素都小于 T[q1] （对应于 x ）， T[q2. .r2] 中的元素都大于等于 T[q1] 。第 11 行计算把输出子数组 A[p3..r3] 分割为 A[p3..q3-1] 和 A[q3+1..r3] 的元素索引 q3 ，接着第 12 行把 T[q1] 直接拷贝到 A[q3] 。

 

       接下来，我们使用嵌套并行进行递归调用。第 13 行 spawn 了对第一个子问题的计算， 14 行并行的进行第 2 个子问题计算。第 15 行的 sync 语句保证过程返回前子问题都已经计算完成。（由于在返回前会隐式执行一条 sync 语句，所以第 15 行的 sync 语句可以省掉，不过明确地写出来是一项好的编码实践）。为了保证在子数组 T[p2..r2] 为空时，代码仍能够正确工作，我们使用了一些巧妙的手段。在每次递归调用中，都会把 T[p1..r1] 的中值元素放到输出子数组中，直到 T[p1..r1] 最后变成空，从而触发基本情形处理。

多线程归并分析

我们先来导出 P-MERGE 的 span P M _∞ (n) ，其中 n 为两个子数组元素的总和： n=n1+n2 。因为第 13 行的 spawn 和第 14 行的调用在逻辑上是并行执行的，所以只需要研究其中开销高一点的那个即可。关键点在于，在最坏的情况下，这两个递归调用中任一个的最多元素个数至多为 3n/4 ，原因如下。因为第 3 到 6 行保证了 n 2 ≤ n1 ，所以有 n2=2n2/2 ≤ (n1+n2)/2 = n/2 。在最坏的情况下，这个递归调用中得有一个要对 T[p1..r1] 中的 └ n1/2 ┘ 个元素和对 T[p2..r2] 中的 全部 n2 元素进行归并，因此该调用中涉及的元素个数为：

 

└ n1/2 ┘+n2 ≤ n1/2 + n2/2 + n2/2

            = ( n1/+ n2)/2 + n2/2

          ≤ n/2 + n/4

          = 3n/4

 

加上第 10 行中 BINARY-SEARCH 调用的开销 Θ(lgn) ，得到最坏情况下， span 的递归式：

 

PM_∞ (n) =  PM_∞ (3n/4) + Θ(lgn)                                             (27.8)

 

( 基础情形的 span 为 Θ(1 ) ，因为第 1 到 8 行的执行时间为常数 ) 。这个递归式不符合 master 定理的任何情形，但是它满足练习 4.6-2 的条件。因此递归式（ 27.8 ）的解为 PM_∞ (n)= Θ(lg² n) 。

 

       现在来分析 P-MERGE 的 work PM₁ (n) ，其值为 Θ(n) 。由于所有 n 个元素都必须得从数组 T 拷贝到数组 A ，因此有 PM₁ (n) = Ω (n) 。接下来，只要证明 PM₁ (n)=O(n) 。

 

       我们首先来导出最坏情况下 work 的递归式。第 10 行的二分查找最坏情况耗时 Θ(lgn) ，其支配了除递归调用之外的其他工作。对于递归调用来说，虽然第 13 行和 14 行的递归调用所归并的元素个数可能不同，但是二者合起来最多归并 n 个元素（实际上是 n-1 个元素，因为 T[q1] 并不参与任何递归调用）。此外，正如我们在分析 span 时所看到的，每个递归调用最多操作 3n/4 个元素。因此得到如下递归式：

 

PM₁ (n) = PM₁ ( ɑ n) + PM₁ ((1- ɑ )n) +O(lgn)                                    (27.9)

 

其中， ɑ位于区间[1/4,3/4] 内，每层递归调用中ɑ的值都有所不同。

 

    我们将通过替换方法来证明递归式（27.9 ）的解为： PM₁ (n)=O(n) 。假设 PM₁ (n ) ≤ c1n-c2lgn （ c1 、 c2 为某个正常量）。通过替换，得到：

 

PM₁ (n ) ≤ ( c1 ɑ n – c2lg( ɑ n) ) +(c1(1- ɑ )n-c2lg((1- ɑ )n))) +Θ(lgn)

      = c1( ɑ +(1- ɑ ) )n-c2(lg( ɑ n)+ lg((1- ɑ) n))+ Θ(lgn)

      = c1n-c2(lg ɑ +lgn+lg(1- ɑ )+lgn)+ Θ(lgn)

      =c1n-c2lgn-(c2(lgn+lg( ɑ (1- ɑ )))- Θ(lgn))

      ≤ c1n-c2lgn

因为我们可以把 c2 选择的足够大，使得 c2(lgn+lg( ɑ (1- ɑ )) ) 支配 Θ(lgn) 项。此外，我们可以把 c1 也选的足够大以满足递归式的基础条件。由于 P-MERGE 的 work PM₁ (n ) 同时为Ω (n) 和 O(n) ，因此 PM₁ (n )= Θ(n) 。

 

       P-MERGE 的 parallelism 为 PM₁ (n )/ PM_∞ (n) = Θ(n/lg² n) 。

多线程归并排序

既然有了一个不错的并行化的多线程归并过程，接下来就可以在多线程归并排序中使用它了。这个版本的归并排序和前面的 MERGE-SORT’ 过程类似，不同之处在于，它多了一个用于存放排序结果的输出子数组 B 作为参数。 P-MERGE-SORT(A,p,r,B,s) 会把 A[p..r] 中的元素进行排序，并把结果保存在 B[s..s+r-p] 中。

 

P-MERGE-SORT(A,p,r,B,s)

1    n =r-p+1

2    if n == 1

3        B[s] = A[p]

4    else 令 T[1..n] 为一个新数组

5        q= └(p+r)/2┘

6        q’=q-p+1

7        spawn P-MERGE-SORT(A,p,q,T,1)

8        P-MERGE-SORT(A,q+1,r,T,q’+1)

9        sync

10       P-MERGE(T,1,q’,q’+1,n,B,s)

 

第 1 行计算输入子数组 A[p..r] 中的元素个数，然后在第 2 、 3 行中处理了当数组仅有一个元素时的基本情形。第 4 到 6 行为第 7 、 8 行中的递归 spawn 和调用（这两者并行运行）做所需要的准备工作。其中，第 4 行创建了一个具有 n 个元素的临时数组 T 来存储递归归并排序的结果。第 5 行计算索引 q ，该索引把 A[p..r] 分割成两个子数组 A[p..q] 和 A[q+1..r] ，后面会递归地对它们进行排序，第 6 行计算第一个子数组 A[p..q] 的元素个数 q’ ，第 8 行使用 q’ 确定 T 中存储排序结果 A[q+1..r] 的起始索引。期间，进行了 spawn 和递归调用，随后是第 9 行的 sync 语句，用来保证被 spawn 的过程执行完毕。最后，第 10 行调用 P-MERGE 对子数组 T[1…q’] 和 T[q’+1..n] 进行归并，把结果放入子数组 B[s..s+r-p] 中。

多线程归并排序分析

首先来分析 P-MERGE-SORG 的 work PMS₁ (n ) ，和 P-MERGE 的 work 比起来要容易分析得多。事实上， work 由如下递归式给出：

 

PMS₁ (n) = 2 PMS₁ (n/2) + PM₁ (n)

       = 2PMS₁ (n/2) + 　 Θ(n)

 

该递归式和 2.3.1 节中普通的 MERGE-SORT 的递归式（ 4.4 ）完全一样，其解为： PMS₁ (n)= Θ(nlgn)( 根据 master 定理中的情形 2) 。

 

       现在来导出并分析最坏情况下的 span PMS_∞ (n) 。因为第 7 、 8 行中的两个对 P-MERGE-SORT 的递归调用在逻辑上是并行的，所以可以忽略其中之一，得到如下递归式：

 

PMS_∞ (n) = PMS_∞ (n/2) + PM_∞ (n)

        = PMS_∞ (n/2) + Θ(lg² n)                                         （ 27.10 ）

 

和递归式（ 27.8 ）一样， master 定理不能应用于递归式（ 27.10 ），但是却符合练习 4.6-2 。其解为 PMS_∞ (n) = Θ(lg³ n) ，所以 P-MERGE-SORT 的 span 为 Θ(lg³ n) 。

 

       并行归并的使用使得 P-MERGE-SORG 在并行度上远远优于 MERGE-SORT’ 。 MERGE-SORT’ 调用了串行 MERGE 过程， parallelism 仅为 Θ(lgn) 。而 P-MERGE-SORT 的 parallelism 为：

 

PMS₁ (n) /PMS_∞ (n) = Θ(nlgn) /Θ(lg³ n) =Θ(nlg² n)

 

无论在理论上还是实践上，这都是一个好得多的结果。在实际实现时，为了降低被渐进表示符号隐藏的那个常量，一种好的做法是通过粗粒度化基础情形来降低一点 parallelsim 。粗粒度化基础情形最直接的方法就是，当数组长度充分小时，就切换到普通的串行排序（比如快速排序）。

练习（略）

问题（略）

本章注解

       并行计算机、并行计算机模型以及并行编程的算法模型已经以不同的形式存在多年了。本书的前面版本中包含了关于排序网络和 PRAM （并行随机访问机器）模型的内容。数据并行模型 [48,168] 是另外一种流行的算法编程模型，向量和矩阵操作的原子化是其主要特征。

 

       Graham[149] 和 Brent[56] 中介绍了 3 个现有的、能够达成定理 27.1 中的界限的调度器， Blumofe 和 Leiserson[52] 中则证明所有贪婪式调度器都能达成这个界限。 Blelloch[47] 提出了一个基于 work 和 span （他称之为计算的“深度”）的算法编程模型，用于数据并行编程。 Blumofe 和 Leiserson[53] 中给出了一个基于随机“工作窃取”的动态多线程的分散式调度算法，并证明它能够达到 E[T _p ] ≤ T ₁ /P+O( T _∞ ) 的边界约束。 Arora ， Blumofe ， Plaxton[19] 以及 Blelloch ， Gibbons ， Matias[49] 中也提出了一些可证明的好算法，用于调度动态多线程计算机。

 

       本章中多线程伪码和编程模型深受 MIT 的 Cilk[51,118] 项目以及由 Cilk Arts 等发布的对于 C++ 的扩展 Cilk++[72] 的影响。本章中的许多多线程算法都来自 C.E.Leiserson 和 H.Prokop 未公开的讲稿，并在 Cilk 或者 Cilk++ 中实现。多线程归并排序算法的灵感来自 Akl[12] 中的一个算法。

 

       顺序一致性的概念来自 Lamport[223] 。