多线程算法（完整版）
——算法导论第3版新增第27章
ThomasH. 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数的递归定义：
F0 = 0,
F1 = 1,
Fi = Fi-1+ Fi-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展示了在计算F6时所创建的递归过程实例树。其中，对于对于FIB(6)的调用会递归地调用FIB(5)和FIB(4)。而对FIB(5)的调用又会去调用FIB(4)。这两个FIB(4)实例返回的结果完全相同（F4=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) ≤ （aFn-1-b）+ (aFn-2-b）+Θ(1)
= a(Fn-1 + Fn-2）- 2b +Θ(1)
= aFn-b – (b-Θ(1))
≤ aFn-b
如果我们在选择b时让其大到足以支配Θ(1)中的常量。那么我们接着可以把a选得大到足以满足初始条件。其分析边界为：
T(n) =Θ(Φn), (27.1)
其中，Φ=(1+sqrt(5))/2是黄金分割率，由等式(3.25)得出。由于Fn随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个处理器上的多线程计算的运行时间。比如，我们用TP来表示算法在P个处理器上的运行时间。work就是在一个处理器上的运行时间，也就是T1。span就是每个strand具有自己独立处理器时的运行时间（也就是说，如果可用的处理器数目是无限的），用T∞来表示。
work和span提供了在P个处理器上运行的多线程计算花费时间TP的下界：
l 在一个单位时间中，具有P个处理器的理想并行计算机最多能够完成P个单位工作，因此在TP时间内，能够完成最多PTP数量的工作。由于总的工作为T1，因此我们有：PTP ≥ T1。两边同除以P得到work法则（work law）：
TP ≥ T1/P. (27.2)
l 具有P个处理器的理想并行计算机肯定无法快过具有无限数量处理器的机器。换种说法，具有无限数量处理器的机器可以通过仅使用P个处理器的方法来仿真具有P个处理器的机器。因此，得到span法则（spaw law）：
TP ≥ T∞. (27.3)
我们用比率T1/ TP来定义在P个处理器上一个计算的加速因子（speedup），它表示该计算在P个处理器上比在1个处理器上快多少倍。根据work法则，TP ≥ T1/P，意味着T1/TP≤P。因此，在P个处理器上的加速因子最多为P。当加速因子和处理器的数目成线性关系时，也就是说，当T1/TP=ΘP时，该计算具有线性加速的性质，当T1/TP=P时，称其为完全的线性加速。
我们把work和span的比率T1/T∞定义为多线程计算的parallelism（并行度）。可以从三个角度来理解parallelism。作为一个比率，parallelism表示了对于关键路径上的每一步，能够并行执行的平均工作量。作为一个上限，parallelism给出了在具有任何数量处理器的机器上，能达到的最大可能加速。最后，也是最重要的，在达成完全线性加速的可能性上，parallelism提供了一个在限制。具体地说，就是一旦处理器的数目超过了parallelism，那么计算就不可能达成完全线性加速。为了说明最后一点，我们假设P > T1/T∞，根据span法则，加速因子满足T1/TP≤T1/T∞<P。此外，如果理想并行计算机的处理器数目P大大超过了parallelism（也就是说，如果P >> T1/T∞），那么T1/TP<<P，这样，加速因子就远小于处理器的数目。换句话说，处理器的数目超过parallelism越多，就越无法达成完全加速。
例如，我们来看看图27.2中P-FIB(4)的计算过程，并假设每个strand花费单位时间。由于work T1=17，span T∞=8，因此parallelism T1/T∞=17/8=2.125。从而，无论我们用多少处理器来执行该计算，都无法获得2倍以上的加速因子。不过，对于更大一些的输入来说，P-FIB(n)会呈现出更大的parallelism。
我们把在一台具有P个处理器的理想并行计算机上执行多线程算法的并行slackness（闲置因子）定义为：(T1/T∞)/P = T1/(PT∞)，也就是计算的parallelism超过机器处理器数目的倍数因子。因此，如果slackness小于1，那么就不能达成完全的线性加速，因为T1/(PT∞)<1，根据span法则，在P个处理器上的加速因子满足T1/TP≤T1/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个处理器上可以达到的最快运行时间为TP= T1/P，根据span法则，最好的情况是TP=T∞。下面的定理表明，因为贪婪式调度器可以以这两个下界之和为其上界，所以其可被证明是一个好的调度器。
定理27.1
在一台具有P个处理器的理想并行计算机上，对于一个wrok为T1，span为T∞的多线程计算，贪婪调度器执行该计算的时间为：
TP≤T1/P + T∞. (27.4)
证明：首先来考虑完全步骤。在每个完全步骤中，P个处理器完成的工作总量为P。我们采用反证法，假设完全步骤的数目严格大于└T1/P┘，那么完全步骤所完成的工作总量至少为：
P*(└T1/P┘+1) = P└T1/P┘ +Ｐ
= T1－（T1 mod P）+ P (根据等式3.8得出)
> T1 (根据不等式3.9得出)。
因此，P个处理器所完成的工作比所需要的还多，矛盾，所以完全步骤的数目最多为└T1/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个处理器的理想并行计算机上，任何由贪婪式调度器调度的多线程计算的运行时间TP，不会超过最优时间的2倍。
证明：令TP*为在具有P个处理器的机器上，一个最优调度器产生的运行时间，令T1和T∞为该计算的work和span。根据work法则和span法则（不等式27.2和27.3），得出：
TP*≥max(T1/P, T∞)，根据定理27.1，有：
TP ≤ T1/P + T∞
≤ 2*max(T1/P, T∞)
≤ 2* TP*
下一个推论告诉我们，对于任何多线程计算来所，随着slackness的增长，贪婪式调度器都可以达到接近完全的线性加速。
推论27.3
令TP为在一台具有P个处理器的理想并行计算机上，贪婪式调度器调度一个多线程计算的运行时间，令T1和T∞为该计算的work和span。那么如果P << T1/T∞，就有TP≈T1/P（或者相等），也就是具有大约为P的加速因子。
证明：假设P<< T1/T∞，那么就有T∞<< T1/P，因此根据定理27.1，有TP≤T1/P + T∞。根据work法则（27.2）得到TP≥T1/P，因此得出TP≈T1/P（或者相等），加速因子为：T1/TP≈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 T1(n)没什么难度，因为我们已经做过了。原始的FIB过程就是P-FIB的串行化版本，因此T1(n)=T(n)= Θ(Φ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为T1 (n)/ T∞ (n) =Θ(Φn/n)，其随着n增长的速度极快。因此，对P-FIB(n)来说，即使在最大的并行计算机上，一个中等大小的n值就足以获得接近完全的线性加速，因为该过程具有相当大的并行slackness。
并行循环
有许多算法，其包含的循环中的所有迭代都是可以并行执行的。我们将看到，可以使用spawn和sync关键字来并行化这种循环，不过如果能够直接指明这种循环的迭代可以并发执行的话，会更加方便一些。我们通过使用parallel并发关键字来在伪码中提供该功能，它位于for循环语句的for关键字之前。
我们以一个n×n的矩阵A=（aij）乘以一个n元向量x=（xj）为例进行说明。相乘的结果为一个n元向量y=（yi），如下：
yi = ∑nj=1aij xj，
i=1，2，…,n。我们可以通过并行地计算y的所有项来进行矩阵-向量的乘法操作，如下：
MAT-VEC(A,x)
1 n = A.rows
2 令y为一个新的长度为n的向量
3 parallel for i = 1 to n
4 yi = 0
5 parallel for i = 1 to n
6 for j = 1 to n
7 yi = yi + aij xj
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 forj = 1 to n
3 yi = yi + aij xj
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 T1(n)，也就是计算其串行化版本的运行时间，这个串行化版本可以通过把parallel for循环替换成普通的for循环得到。由此，我们得到T1(n)= Θ(n2)，因为第5到7行的两重嵌套循环所产生的平方级运行时间占支配地位。在这个分析中，我们忽略掉了实现并行循环的递归spawn的开销。事实上，和其串行化版本相比，递归spawn的开销确实增加了并行循环的工作量，不过并不是渐进关系的。原因如下，因为递归过程实例树是一颗满二叉树，所以内部节点的个数正好比叶子的个数少1（见练习B.5-3）。每个内部节点分割迭代范围时所耗费的都是常数时间，并且每个叶子都对应循环中的一个迭代，其至少耗费常数时间（在本例中是Θ(n)）。因此，我们可以把递归spawn的开销分摊到迭代的工作中，对全部工作来说，至多增加了一个常数倍数因此。
在实际实现中，动态多线程并发平台时常会在一个叶子中执行多个迭代，从而使得递归产生的叶子的粒度变粗，这个过程可以是自动地，也可以由程序员来控制，因此减少了递归spawn的开销。付出的代价是降低了并行度，不过，如果计算具有局够大的并行slackness，那么还是可以达成接近完全的线性加速的。
在分析并行循环的span时，也必须得考虑到递归spawn的开销。由于递归调用的深度和迭代的次数成对数关系，因此对于一个具有n次迭代，其第i个迭代的span为iter∞ (i)的并行循环来说，其span为：
T∞(n)= Θ(lgn)+ max1≤i≤niter∞ (i)。
例如，对于以一个n×n矩阵为参数的MAT-VEC来说，第3、4行中的并行初始化循环的span为Θ(lgn)，因为和每个迭代中的常数工作时间相比，递归spawn占支配地位。第5到7行中的双重嵌套循环的span为Θ(n)，因为外层parallel for循环的每个迭代都包含着内层（串行）for循环的n个迭代。伪码中剩余部分的span为常数，因此整个过程的span由双重嵌套循环支配，也就是Θ(n)。由于过程的work为Θ(n2)，所以parallelism为Θ(n2)/ Θ(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 yi = 0
5 parallel for i =1 to n
6 parallelfor j = 1 to n
7 yi = yi + aij xj
8 return y
遗憾的是，这个过程是错误的，因为在第7行更新yi时存在竞争问题（对于j的n个值并行去更新）。练习27.1-6会让读者提供一个具有Θ(lgn)的span值的正确实现。
具有竞争问题的多线程算法有时也是正确的。例如，两个并行线程可以把相同的值存入一个共享的变量中，谁先谁后无关紧要。不过我们通常认为具有竞争问题的代码都是非法的。
来自国际象棋程序的教训
我们以一个真实的故事来结束本小节，该故事发生在世界级多线程国际象棋对弈程序Socrates[81]的开发期间（时间做了简略）。该程序的原型是在一个具有32个处理器的计算机上进行的，但是最终运行在一个具有512个处理器的超级计算机上。某天，开发人员对程序进行了一项优化，针对一项重要的在32个处理器上基准测试（benchmark），把其运行时间从 T32= 65秒减少到T’32=40秒。然而，开发人员使用关于work和span的性能度量得出结论：在32个处理器上更快的优化版本，在512个处理器上运行时，实际上比原始版本慢一些。结果，他们放弃了“优化”。
他们的分析方法是这样的。程序的原始版本的work T1= 2048秒，span T∞=1秒。如果把不等式(27.4)看做是等式，TP=T1 /P+T∞，并把它当作在P个处理器上的近似运行时间，确实可以得出，T32=2048/32+1＝ 65。如果做了优化，work变成T’1=1024秒，span变成T’∞= 8秒。采用同样的近似方法，有T’32=1024/32＋8=40。
但是，当在512个处理器上进行计算运行时间时，这两个版本的相对速度会发生转换。此时，T512=2048/512+1＝5秒，T’512=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 parallelfor i = 1 to n
4 parallelfor j = 1 to n
5 cij = 0
6 for k = 1 to n
7 cij = cij+ aik* bkj
8 returnC
现在来分析这个算法，由于该算法的串行化版本就是SQUARE-MATRIX-MULTIPLY，因此其work为T1(n)= Θ(n3)，和SQUARE-MATRIX-MULTIPLY的运行时间一致。其span为T∞(n)= Θ(n)，因为执行过程先沿着第3行启动的parallel for循环所产生的递归树向下，然后再沿着第4行启动的parallel for循环所属产生的递归树向下，接着执行了第6行中普通for循环中的所有n个迭代，从而总的span为，Θ(lgn)+ Θ(lgn)＋Θ(n)＝Θ(n)。因此，其parallelism为Θ(n3)/ Θ(n)=Θ(n2) 。练习27.2-3会要求读者并行化内层循环得到一个span为Θ(lgn)的算法，不能直接使用parallel for，因为会产生竞争问题。
矩阵相乘的多线程分治算法
如4.2节中所讲，可以采用Strassen分治策略在Θ(nlg7)= Θ(n2.81)的时间内串行地完成n×n矩阵的相乘，这激发我们去寻找该算法的多线程版本。和4.2节中一样，我们先来多线程化一个简单一些的分治算法。
回忆一下第77页中的SQUARE-MATRIX-MULTIPLY-RECURSIVE过程，它把两个n×n矩阵A、B相乘得到一个n×n矩阵C，采用的方法是把这三个矩阵都分割成四个n/2×n/2的子矩阵：
A={{A11,A12},{ A21,A22}},B={{B11,B12},{ B21,B22}}, C={{C11,C12},{C21,C22}}。
我们可以把矩阵积记为：
{{C11,C12},{ C21,C22}} = {{A11,A12},{A21,A22}} * {{B11,B12},{ B21,B22}}
= {{A11B11, A11 B12},{ A21 B11,A21 B12}} +
{{A12B21, A12 B22},{ A22 B21,A22 B22}} (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 ifn == 1
3 c11 = a11 b11
4 else 令T为新的n×n矩阵
5 把A、B、C和T分割成n/2×n/2的子矩阵
A11、A12、A21、A22；B11、B12、B21、B22；C11、C12、C21、C22；
以及T11、T12、T21、T22
6 spawn P-MATRIX-MULTIPLY-RECURSIVE(C11, A11, B11)
7 spawn P-MATRIX-MULTIPLY-RECURSIVE(C12, A11, B12)
8 spawn P-MATRIX-MULTIPLY-RECURSIVE(C21, A21, B11)
9 spawn P-MATRIX-MULTIPLY-RECURSIVE(C22, A21, B12)
10 spawn P-MATRIX-MULTIPLY-RECURSIVE(T11, A12, B21)
11 spawn P-MATRIX-MULTIPLY-RECURSIVE(T12, A12, B22)
12 spawn P-MATRIX-MULTIPLY-RECURSIVE(T21, A22, B21)
13 P-MATRIX-MULTIPLY-RECURSIVE(T22, A22, B22)
14 sync
15 parallel for i = 1 to n
16 parallel for j = 1 to n
17 cij= cij+ tij
第3行对基本情形进行了处理，也就是进行1×1矩阵的相乘。第4到17行处理了递归情形。在第4行中创建了一个临时矩阵T，在第5行中吧矩阵A、B、C和T分割成n/2×n/2的子矩阵。（和第77页的SQUARE-MATRIX-MULTIPLY-RECURSIVE一样，我们忽略如何用索引计算来表示矩阵的子矩阵部分这个小问题）。第6行的递归调用把子矩阵C11设置为子矩阵A11 和B11乘积，这样C11就等于等式（27.6）中所形成其和的两项中的第一个。同样，第7到9行把C12、C21和C22设置为等于等式（27.6）中形成其和的两项中的第一项。第10行把子矩阵T11设置为子矩阵A12和B21乘积，这样T11就等于形成C11和的两项中的第二个。第11到13行分别把T12、T21和T22设置为形成C12、C21、C22和的两项中的第二个。前7个跌鬼调用时spawn出来的，最后一个运行在主strand中。第14行中sync语句保证了第6到13行中的矩阵积的计算都已经完成，然后在第15到17行中，使用双重嵌套parallelfor循环把积T和C相加。
首先，仿照其原始SQUARE-MATRIX-MULTIPLY-RECURSIVE过程的串行化运行时间分析方法，来分析一下P-MATRIX-MULTIPLY-RECURSIVE过程的work M1(n)。在递归情形中，分割用时为Θ(1)，然后执行8个n/2×n/2矩阵的递归乘法，最后耗时Θ(n2)进行两个n×n矩阵相加。因此，其work M1(n)的递归等式为：
M1(n) =8 M1(n/2) +Θ(n2)
=Θ(n3) （根据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)= Θ(lg2n)。
既然知道了P-MATRIX-MULTIPLY-RECURSIVE的work和span，就可以计算其parallelism，为M1(n)/ M∞(n) =Θ(n3/lg2n)，非常之高。
Strassen方法的多线程化
我们根据和第79页中相同的方法来多线程化Strassen算法，仅使用嵌套并行：
1、 把输入矩阵A、B以及输出矩阵按照等式（27.6）分割成n/2×n/2的子矩阵。使用索引计算，该步骤的wiork和span均为Θ(1)。
2、 创建10个矩阵S1、S2、…、S10，每个都是n/2×n/2的，值为第一步中创建矩阵的和或者差。使用双重嵌套parallel for循环，创建这10个矩阵的work和span分别为：Θ(n2)和Θ(lgn)。
3、 使用第1步中创建的子矩阵和第2步中创建的10个矩阵，递归地spawn出7个计算来计算7个n/2×n/2的矩阵积P1、P2、…、P7。
4、 通过对Pi矩阵不同组合进行加、减，同样采用双重嵌套parallel loop循环计算出所期望的结果矩阵C的子矩阵C11、C12、C21、C22。计算着四个子矩阵的work和span分别为：Θ(n2)和Θ(lgn)。
现在来分析这个算法，由于其串行化版本和原始的串行算法一样，因此work就是串行化的运行时间，也就是Θ(nlg7)。对于P-MATRIX-MULTIPLY-RECURSIVE来说，我们可以设计出关于span的递归等式。在本例中，虽然7个递归调用并行执行，但是由于它们操作的矩阵大小完全一样，从而可以像P-MATRIX-MULTIPLY-RECURSIVE中一样得到相同的递归式(27.7)，其解为Θ(lg2n)。因此，多线程Strassen算法的parallelsism为Θ(nlg7)/ Θ(lg2n)，虽然比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’1(n)：
MS’1(n)= 2MS’1(n/2) +Θ(n) = Θ(nlgn)，
和归并排序的串行执行时间一样。由于对MERGE-SORT’的两个递归调用可以并行运行，因此span MS∞’(n)由如下递归式定义：
MS∞’(n) = MS∞’(n/2)+Θ(n) ＝ Θ(n)。
因此，MERGE-SORT’的parallelism为MS’1(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]一样长，也就是说，n1≥n2，所以在基本情形中只需检查是否n1=0即可。我们必须得保证，在递归中正确地处理了两个子数组中只有一个为空的情况，根据我们n1≥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行确保n1≥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 PM∞(n)，其中n为两个子数组元素的总和：n=n1+n2。因为第13行的spawn和第14行的调用在逻辑上是并行执行的，所以只需要研究其中开销高一点的那个即可。关键点在于，在最坏的情况下，这两个递归调用中任一个的最多元素个数至多为3n/4，原因如下。因为第3到6行保证了n2≤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)= Θ(lg2n)。
现在来分析P-MERGE的work PM1(n)，其值为Θ(n)。由于所有n个元素都必须得从数组T拷贝到数组A，因此有PM1(n) =Ω(n)。接下来，只要证明PM1(n)=O(n)。
我们首先来导出最坏情况下work的递归式。第10行的二分查找最坏情况耗时Θ(lgn)，其支配了除递归调用之外的其他工作。对于递归调用来说，虽然第13行和14行的递归调用所归并的元素个数可能不同，但是二者合起来最多归并n个元素（实际上是n-1个元素，因为T[q1]并不参与任何递归调用）。此外，正如我们在分析span时所看到的，每个递归调用最多操作3n/4个元素。因此得到如下递归式：
PM1(n) = PM1(ɑn) + PM1((1-ɑ)n) +O(lgn) (27.9)
其中，ɑ位于区间[1/4,3/4]内，每层递归调用中ɑ的值都有所不同。
我们将通过替换方法来证明递归式（27.9）的解为：PM1(n)=O(n)。假设PM1(n)≤c1n-c2lgn（c1、c2为某个正常量）。通过替换，得到：
PM1(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 PM1(n)同时为Ω(n)和O(n)，因此PM1(n)= Θ(n)。
P-MERGE的parallelism为PM1(n)/PM∞(n) =Θ(n/lg2n)。
多线程归并排序
既然有了一个不错的并行化的多线程归并过程，接下来就可以在多线程归并排序中使用它了。这个版本的归并排序和前面的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 PMS1(n)，和P-MERGE的work比起来要容易分析得多。事实上，work由如下递归式给出：
PMS1(n) = 2 PMS1(n/2) + PM1(n)
= 2PMS1(n/2) + Θ(n)
该递归式和2.3.1节中普通的MERGE-SORT的递归式（4.4）完全一样，其解为：PMS1(n)= Θ(nlgn)(根据master定理中的情形2)。
现在来导出并分析最坏情况下的span PMS∞(n)。因为第7、8行中的两个对P-MERGE-SORT的递归调用在逻辑上是并行的，所以可以忽略其中之一，得到如下递归式：
PMS∞(n)= PMS∞(n/2) + PM∞(n)
= PMS∞(n/2) +Θ(lg2n) （27.10）
和递归式（27.8）一样，master定理不能应用于递归式（27.10），但是却符合练习4.6-2。其解为PMS∞(n) =Θ(lg3n)，所以P-MERGE-SORT的span为Θ(lg3n)。
并行归并的使用使得P-MERGE-SORG在并行度上远远优于MERGE-SORT’。 MERGE-SORT’调用了串行MERGE过程，parallelism仅为Θ(lgn)。而P-MERGE-SORT的parallelism为：
PMS1(n)/PMS∞(n) =Θ(nlgn) /Θ(lg3n) =Θ(nlg2n)
无论在理论上还是实践上，这都是一个好得多的结果。在实际实现时，为了降低被渐进表示符号隐藏的那个常量，一种好的做法是通过粗粒度化基础情形来降低一点parallelsim。粗粒度化基础情形最直接的方法就是，当数组长度充分小时，就切换到普通的串行排序（比如快速排序）。
练习（略）
问题（略）
本章注解
并行计算机、并行计算机模型以及并行编程的算法模型已经以不同的形式存在多年了。本书的前面版本中包含了关于排序网络和PRAM（并行随机访问机器）模型的内容。数据并行模型[48,168]是另外一种流行的算法编程模型，向量和矩阵操作的原子化是其主要特征。
Graham[149]和Brent[56]中介绍了3个现有的、能够达成定理27.1中的界限的调度器，Blumofe和Leiserson[52]中则证明所有贪婪式调度器都能达成这个界限。Blelloch[47]提出了一个基于work和span（他称之为计算的“深度”）的算法编程模型，用于数据并行编程。Blumofe和Leiserson[53]中给出了一个基于随机“工作窃取”的动态多线程的分散式调度算法，并证明它能够达到E[Tp]≤T1/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]。
