并行程序设计导论

1 为什么要并行化

1.1 为什么需要不断提升的性能

• 不断提升的计算能力称为许多飞速发展领域的核心力量
• 气候模拟
• 蛋白质折叠
• 药物发现
• 能源研究
• 数据分析

1.2 为什么需要构建并行系统

• 单处理器性能提升的主要原因之一：日益增加的集成电路晶体管密度。
• 通过继续增快继承电路的速度来提高处理器性能的方法变得不再可行。
• 如何利用不断增加的晶体管密度：并行。
• 集成电路制造商：在单个芯片上放置多个处理器。（多核处理器）

1.3 为什么需要编写并行程序

• 串行程序改写为并行程序。
• 最好的并行化，可能通过一步步回溯，然后发现一个全新的算法来得到。

代码分析：计算 n 数的值再累加求和，串行代码：

sum = 0;
for(i=0; i<n; i++){
    x = Compute_next_value(...);
    sum += x;
}

假设有 p 个核，且 p 远小于 n，那么每个核能够计算大约 n/p 个数的值并累加求和，以得到部分和：

my_sum = 0;
my_first_i = ...;
my_last_i = ...;
for(my_i = my_first_i; my_i < my_last_i; my_i++){
    my_x = Compute_next_value(...);
    my_sum += my_sum;
}

此块的前缀 my_ 代表每个核都使用自己的私有变量，并且每个核能够独立于其他核来执行该代码块。

方法1 每个核计算完各自的 my_sum 值后，将自己的结果值发送给一个指定为 “master” 的核（主核），master 核收到的部分和累加得到全局和：

if(I'm the master core){
    sum = my_sum;
    for each core other than myself{
        receive value from core;
        sum += value;
    }
}else{
    send my_sum to the master;
}

方法2 当核的数目比较多时，不再由 master 核计算所有部分和的累加工作，可以两两结对，0核将自己的部分和与 1 号核的部分和做加法，2号核将自己的部分和与3号核的部分和做加法，4号核将自己的部分和与5号核的部分和做加法，以此类推。然后，再在偶数核上重复累加部分和：0号核加上2号核，4号核加上6号核，以此类推。

如图1-1所示。圆圈表示当前核所得到的和，箭头表示一个核将自己的部分和发送给另一个核，加号表示一个核在收到另一个核发送来的部分和后与自己本身的部分相加。

分析：上述两种计算全局总和的算法中，master 核承担了更多的工作量，整个程序计算全局总和的时间 = master 核的计算时间。 在 8 个核的情况下，① 第一种方法，master 核需要执行 7 次接收操作。② 第二种方法，master 核仅需执行 3 次执行操作，提高了 2 倍。 在 1000 个核的情况下，① 第一种方法需要 999 次接收和加法操作，② 第二种方法 10次，提高了100倍。

1. 第一种计算全局总和的方法是对串行求和程序的一般化。
2. 第二种计算全局总和的方法与原来的串行程序并没有多大关系。

1.4 怎样编写并行程序

大部分方案的基本思想：将要完成的任务分配给各个核。

两种广泛采用的方法：

• 任务并行：将待解决问题所需要执行的各个任务分配到各个核上执行。
• 数据并行：将待解决问题所需要处理的数据分配给各个核，每个核在分配到的数据集上执行大致相似的操作。

计算全局总和的例子：

• 第一部分可以认为是 数据并行 的一个实例。数据是通过 Compute_next_value 计算得到的值，该值在所赋予的数据集上执行大致相同的操作。
• 第二部分认为是 任务并行 的一个实例。总共有两个任务：
  1. 一个任务由 master 核执行，负责接收 其他核 传来的 部分和，并 累加部分和。
  2. 另一个任务由 其他核 执行，负责将自己计算得到的部分和传递给 master 核。
• 协调过程 包括：
  • 通信，一个或多个核将自己的 部分和结果 发送给其他的核。
  • 负载均衡，在各个核心之间均匀分配工作，以避免某个核心过载。
  • 同步，由于每个核心以自己的速度工作，确保各核心之间不会相差太远。

1.5 我们将做什么

两种并行系统：

• 共享内存系统
  • 每个核能够共享访问计算的内存，理论上每个核能够读、写内存的所有区域。因此可以通过检测和更新共享内存的数据来协调各个核。
• 分布式内存系统
  • 每个核都拥有自己的私有内存，核之间的通信是显式的，必须使用类似于在网络中发送消息的机制。

1.6 并发、并行、分布式

• 并发计算：一个程序的多个任务在同一时间段内可以同时执行。
• 并行计算：一个程序通过多个任务紧密协作来解决某个问题。
• 分布式计算：一个程序与其他程序协作来解决某个问题。

并行程序和分布式程序都是并发的。

2 并行硬件和并行软件

2.1 背景知识

冯·诺依曼结构

• 冯·诺依曼结构：主存、中央处理单元、以及 主存 和 CPU 之间的互连结构。
  • 主存：有很多区域，每个区域都可以存储指令和数据。每个区域都有一个地址，可以通过这个地址来访问向相应的区域及区域中存储的数据和指令。
  • 中央处理单元 (CPU)：分为 控制单元 和 算术逻辑单元。
    • 控制单元：决定应该执行程序中的哪些指令。
    • 算术逻辑单元：负责执行指令。
    • 寄存器：存储 CPU 中的数据和程序执行时的状态信息。
    • 程序计数器：用于存放下一条指令地址的特殊的寄存器。
• 数据/指令传送到 CPU = 数据/指令从内存中取出或读出。
• 输出/指令从CPU 传送到 主存中 = 数据/指令 写入或存入 内存中。
• 冯·诺依曼瓶颈：主存和 CPU 之间的分离带来的数据访问瓶颈。

进程、多任务及线程

• 操作系统：用来管理计算机的软件和硬件资源的主要软件。

• 进程：运行着的程序的一个实例。

  一个进程包括如下实体：

  • 可执行的机器语言程序。
  • 一块内存空间。
  • 操作系统分配给进程的资源描述符。
  • 安全信息。
  • 进程状态信息。

2.2 对冯诺依曼模型的改进

• Cache (缓存)

  高速缓冲存储器 (cache, 简称 缓存) 的访问时间 比其他存储区域的访问时间短。

  CPU 缓存 是一组相比于主存，CPU 能更快访问的 内存区域。

  • 空间局部性：程序访问完一个存储区域往往会访问接下来的区域。

  • 局部性：在访问完一个内存区域（指令或者数据），程序会在不久的将来（时间局部性）访问邻近的区域（空间局部性）。

  • 为了运用局部性原理，系统使用 更宽的互连结构 来访问数据的指令。

  • 一次内存访问能存取一整块代码和数据，而不只是 单条指令 和 单条数据。这些块称为 高速缓存块 或 高速缓存行。

  • Cache 分为不同的层。

  • 如果 CPU 需要访问指令或数据时，它会沿着 Cache 的层次结构向下查询，首先查询 L1 cache，接着 L2 cache。以此类推，如果 Cache 中没有所需要的信息，就会访问主存。

  • 当向 Cache 查询信息时，Cache 中有信息，则称为 Cache命中 或 命中。

    如果信息不存在，则称为 Cache缺失 或 缺失。

  • 当 CPU 向 Cache 中写数据时，Cache 中的值就会 不同 或 不一致。解决不一致性的方法：

    • 写直达 Cache中：当 CPU 向 Cache 写数据时，高速缓存会立即存入主存。
    • 写回 Cache 中：数据不是立即更新到主存中，而是将发生数据更新的高速缓存行标记成脏，当发生高速缓存替换时，标记为脏的高速缓存行倍写入主存。

• Cache 映射

  • 全相联 Cache：每个高速缓存行能够放置在 Cache 中的任意位置。
  • 直接映射 Cache：每个高速缓存行在 Cache 中有唯一的位置。
  • n 路组相联：每个高速缓存能放置到 Cache 中 n 个不同区域位置的一个。
  • 当内存中的行（多于一行）能被映射到 Cache 中的多个不同位置时，需要决定 替换 或者 驱逐 Cache 中的哪一行。
  • 最常用的替换方案：最近最少使用，CPU 记录各个块的被访问的次数，替换最近访问次数最少的块。

• Cache 和程序：一个实例

  C 语言以 “行主序” 来存储二维数据，行主序存储模式下，先存储二维数组的第一行，紧接着第二行，以此类推。

  

  分析：上面代码中，第一个嵌套循环 比 第二个嵌套循环 具有更好的性能，因为它 顺序访问二维数组 中的数据。

• 虚拟存储器

  • 利用 虚拟存储器，主存可以作为辅存的缓存。它通过在主存中存放当前执行程序所需要用的部分，利用时间和空间局部性；那些暂时用不到的部分存储在辅存的块中，称为 交换空间。
  • 虚拟存储器也是对数据块和指令块进行 操作。这些块通常称为页。
    • 页，通常很大。
    • 大多数系统采用固定大小的页，从 4~16kB 不等
  • 虚拟页号：
    • 当程序被编译时，其页面会被分配虚拟页号。
    • 当程序运行时，会创建一个表，将虚拟页号映射到物理地址。
    • 使用页表将虚拟地址转换为物理地址。
  • 使用页表会增加程序总体的运行时间。
  • 处理器有一种专门用于地址转换的缓存：转译后备缓冲区 (Translation-Lookaside Buffer, TLB)：在快速存储介质中缓存一些页表的条目。
    • 如果想要访问的页面不在内存中，即页表中该页没有合法的物理地址，该页只存储在磁盘上，那么这次访问称为 页面失效。
    • 磁盘访问代价太大，应该尽可能避免，所以虚拟内存通常采用写回方案。
    • 采用 写回 方案。

• 指令集并行

  • 指令集并行：通过让多个处理器部件或者功能单元同时执行指令来提高处理器的性能。

  • 两种主要方法实现指令集并行：

    • 流水线：将功能单元分阶段安排。
    • 多发射：让多条指令同时启动。

  • 流水线

    

    • 每个操作 1 ns，则加法操作 7 ns，这个 for 循环需要 7000 纳秒。
    • 另一个方案：将浮点数加法器划分为 7 个独立的硬件或功能单元。第一个单元取两个操作数，第二个比较指数，以此类推。假设一个功能单元的输出是下面一个功能单元的输入，那么加法功能单元的输出是规格化结果的功能单元的输入。在时间5后，流水循环每1纳秒产生一个结果，而不再是每 7 纳秒一次，所以，执行 for 循环的总时间从 7000 纳秒降低到 1006 纳秒，提高了近 7 倍。

  • 多发射

    • 流水线通过将功能分成多个单独的硬件或功能单元，并把它们按顺序串接起来提高性能。
    • 多发射处理器通过 复制功能单元 来同时执行程序中的不同指令。
    • 静态多发射：功能单元编译时调度
    • 动态多发射：运行时间调度的。
    • 超标量：一个支持 动态多发射的处理器。
    • 为了能够利用 多发射，系统必须找出能够同时执行的指令，其中一种最重要的技术是 预测。
    • 在预测技术中，编译器或者处理器对一条指令进行猜测，然后再猜测的基础上执行代码。

  • 硬件多线程

    • 线程级并行 (Thread-Level Parallelism, TLP)：尝试通过同时执行不同线程来提供并行性。TLP 提供的是粗粒度并行，即同时执行的程序基本单元（线程）比细粒度的程序单元更大。

    • 硬件多线程：为系统提供了一种机制，使得当前执行的任务被阻塞，系统能够继续其他有用的工作。

    • 细粒度 多线程中，线程在每条指令执行完后切换线程，而跳过被阻塞的线程。

    • 同步多线程：通过允许 多个线程同时使用多个功能单元 来利用超标量处理器的性能。

      我们指定“优先”线程，那么能够在一定程度上减轻线程减速的问题。优先线程 是指有多条指令就绪的线程。

2.3 并行硬件

2.3.1 SIMD 系统

• Flynn 分类法：用来对 计算机体系结构 进行分类。按照它能够同时管理的指令流数目和数据流数目来对系统分类。

• 单指令多数据流 (SIMD)系统是 并行系统。

  • 将数据项分配给各个处理器来实现并行性。（数据并行）
  • SIMD 系统通过对多个数据执行相同的指令从而实现在多个数据流上的操作。
  • 缺点：所有算术逻辑单元必须执行相同的指令，否则保持空闲。在经典的设计中，它们还必须同步操作。ALU 没有指令存储，对于大规模数据并行问题非常有效，但对于其他类型更复杂的并行问题则不然。

• 向量处理器：能够对数组或者数据向量进行操作。

  典型的特性有以下特征：

  • 向量寄存器。能够存储由多个操作数组成的向量，同时能够同时对其内容进行操作的寄存器。
  • 向量化和流水化的功能单元。对向量中的每个元组做同样的操作。
  • 向量指令。在向量上操作的指令。
  • 交叉存储器：由多个内存“体”组成，每个内存体能够独立访问。在访问完一个内存体之后，再次访问它之前需要一个时间延迟，但如果接下来的内存访问是另一个内存体，那么它很快能访问到。
  • 步长式存储器访问和硬件散射/聚集：程序能够访问向量中固定间隔的元组。
  • 速度快，容易使用，高存储带宽。
  • 不能处理不规则的数据结构和其他的并行结构。可扩展性是个限制。可扩展性：能够处理更大问题的能力。

• 图形处理单元 (GPU)

  • 实时图形应用编程接口使用 点、线、三角形来表示物体的表面。

  • 使用图形处理流水线 把物体表面的内部表示 转换为一个像素的数组。

  • 编程阶段的行为可以通过 着色函数 来说明。

  • GPU 通过使用 SIMD 并行来优化性能。

    这是通过在每个 GPU 处理核引入大量的 ALU 来获取的。

  • 实际线程的数目依赖于着色函数需要的资源的数量。

2.3.2 MIMD 系统

• 系统 支持同时多个指令流在多个数据流上操作。

  通常包括一组完全独立的处理单元或者核，每个处理单元或者核都有自己的控制单元和 ALU。

• MIMD 通常是异步的。

• MIMD 主要有两种类型：共享内存系统 和 分布式内存系统。

• 共享内存系统：一组自治的处理器通过互联网络与内存系统相互连接，每个处理器能够访问每个内存区域。在共享内存系统中，处理器通过访问共享的数据结构来隐式地通信。

  • 广泛使用的共享内存系统使用一个或多个多核处理器。一个多核处理器在一块芯片上有多个 CPU 或核。通常，每个核都拥有私有的 L1 Cache，而其他的 Cache 可以在核之间共享，也可以不共享。

  • 在拥有多个多核处理器的共享内存系统中，共享内存的方式：

    • 第一种：互联网络将所有的处理器连接到主存。每个核访问内存中任何一个区域的时间都相同。“一致内存访问 (UMA)”

      容易编程。

    • 第二种：将每个处理器连接到一块内存，通过处理器中内置的特殊硬件使得各个处理器可以访问内存中的其他快。访问与核直接连接的那块内存区域比访问其他内存区域要快很多。“非一致内存访问 (NUMA)”

      可以使用更大的内存。

• 分布式内存系统：每个处理器都有自己私有的内存空间，处理器-内存对之间通过互联网络相互通信。处理器之间是通过发送消息或者使用特殊的函数来访问其他处理器的内存，从而进行显式的通信。

  • 最广泛使用的分布式内存系统：集群 (clusters)，它们由一组商品化系统组成，通过商品化网络连接。实际上，这些系统中的节点（通过通信网络相互连接的独立计算单元），通常都是一个或多个多核处理器的共享内存系统。这种系统通常称为混合系统。
  • 网络 提供一种基础架构，使地理上分布的计算机大型网络转换成一个分布式内存系统。通常，这样的系统是异构的，即每个节点都是由不同的硬件构造的。

2.3.4 Cache 一致性

• 情景分析

  假如共享系统中有两个核，每个核都有各自的私有的数据 Cache。

  假设，x 是一个共享变量并初始化为 2，y0 是 核0 私有的，y1 和 z1 是 核1 私有的。现在假设按指定时序执行下面的语句：

  时间	核0	核1
  0	y0 = x;	y1 = 3*x;
  1	x = 7;	没有包含x的语句
  2	没有包含x的语句	z1 = 4*x

  • y0 的内存区域最终得到值 2，y1 的内存区域最终得到值6。
  • z1 的值无法确定：时间 0 时刻，x 已经在核 1 的 Cache 中。
    • 除非出于某些原因，x 清除出核0的 Cache，然后再重新加载到核1 的 Cache 中，此时 z1 赋值使用的是最新值 x = 7，z1 = 28。
    • 否则，z1 仍然使用 Cache 中的 x = 2，得到 z1 = 8。

• 这种不可预测的行为与系统使用写直达还是写回策略无关。

  核0 在 时间1 更新主存的时候，并不会修改 核1 的 Cache 值。

• 两种主要的方法来保证 Cache 的一致性：监听 Cache 一致性协议 和 基于目录的 Cache 一致性协议。

• 监听 Cache 一致性协议

  当一个核更新Cache 中 x 的副本，就会广播通知其他核包含 x 的整个 Cache 行已经更新。

  1. 互联网络不一定是总线，只要能够支持从每个处理器广播到其他处理器。
  2. 如果是写直达 Cache，每个核都能“检测”写。
  3. 如果是写回 Cache，就需要额外的通信，因为对 Cache 的更新不会立即发送给内存。

• 基于目录的 Cache 一致性协议

  使用一个叫 目录 (directory) 的数据结构。目录存储每个内存行的状态。一般地，这个数据结构是分布式的。

  • 假设 每个核/内存对 负责存储一部分的目录。这部分目录表示局部内存对应高速缓存行的状态。
  • 当一个高速缓存行被读入时，与这个高速缓存行相对应的目录项就会更新，表示 核 0 有这个行的副本。当一个变量需要更新时，就会查询目录，并将所有包含该变量高速缓存行置为非法。
  • 目录需要大量额外的存储空间。

• 伪共享

  • 原因：CPU Cache 是由硬件来实现的。因为 硬件 是对 高速缓存行 进行操作而不是对单独的变量进行操作。这个特点可能会对性能带来极坏的影响。

  • 代码分析

    C

    折叠 复制

    int i, j, m, n;
    double y[m];
    /* 赋值 y = 0 */
    ...
    for(i = 0; i < m; i++)
        for(j = 0; j < n; j++)
            y[i] += f(i, j);

    重复调用 f(i, j) ，并将计算的值添加到一个向量中。

    为了并行化，将外部循环的各次迭代分配给各个处理器核来处理。假如有 core_count 个核，可能把第一个 m/core_count 个迭代分配给第一个核，下一个 m/core_count 个迭代分配给第二个核，以此类推。

    C

    折叠 复制

    /* 私有变量 */
    int i, j, iter_count;
    /* 共享变量 */
    int m, n, core_count;
    double y[m];
    
    iter_count = m/core_count;
    
    /* 核0 进行以下操作 */
    for(i = 0; i < iter_count; i++)
        for(j = 0; j < n; j++)
            y[i] += f(i, j);
    
    /* 核1 进行以下操作 */
    for(i = iter_count + 1; i < 2*iter_count; i++)
        for(j = 0; j < n; j++)
            y[i] += f(i, j);
    ...

    现在假设 共享内存系统 有 两个核，m=8 ，双精度浮点数长度为 8 ，高速缓存行是 64字节，y[0] 存储在高速缓存行的起始位置。一个高速缓存行能够存储 8 个双精度浮点数，y 占一整个高速缓存行中。

    • 当核 0 和 核1 同时执行代码时，因为 y 的所有值都存储在单个高速缓存行中，所以每次一个核执行语句 y[i] += f(i, j) 时，高速缓存行就会失效，
    • 当另一个核尝试执行该语句时，就必须从内存中将更新过的高速缓存行取出。所以，如果 n 很大，尽管 核0 和 核1 不会访问 y 中对方的元素，但我们能预料到，大量的赋值操作 y[i] += f(i, j) 会访问主存。这称为 伪共享。因为系统表现的好像核之间会共享 y。

    多个线程 同时修改位于 同一缓存行 中的 不同变量，导致缓存行频繁失效，进而降低性能。

    伪共享 不会引发错误结果，但是，它能引起过多不必要的访存，降低程序的性能。可以通过在线程或者进行中临时存储数据，再把临时存储的数据更新到共享存储来降低共享带来的影响。

2.3.5 共享内存与分布式内存

• 分布式内存系统比较适合于那些需要大量数据和计算的问题。

2.4 并行软件

在运行共享内存系统中，会启动一个进程，然后派生 (folk) 出多个线程。

讨论共享内存程序时，通常指 正在执行任务的线程。

当讨论分布式内存程序时，使用的是多个处理器，所以我们指的是正在执行任务的进程。

2.4.1 注意事项

本书余下部分只讨论 MIMD 系统的软件。

• 主要关注 单程序多数据流 (Single Program, Multiple Data, SPMD) 程序。

• SPMD 程序 仅包含一段可执行代码，通过条件转移语句，可以让一段代码在执行表现得像是在不同处理器上执行不同的程序。

  PLAINTEXT

  折叠 复制

  if(I'm thread/process 0)
      do this;
  else
      do that;

• SPMD 程序也能够实现数据并行：

  PLAINTEXT

  折叠 复制

  if(I'm thread/process 0)
      operate on the first half of the array;
  else
      operate on the second half of the array;

• 任务并行：将任务划分，分给各个进程或线程来实现并行。

2.4.2 进程或线程的协调

假如有两个数组，我们想要将它们相加：

double x[n], y[n];

for(int i = 0; i < n; i++)
    x[i] += y[i];

为了并行化这段代码，只需要将数组中的元素分配给线程或者进程。

在这个例子，程序只需要

1. 将任务在 进程/线程 之间分配

   ① 这个分配可以使得每个进程/线程 获得 大致相等的工作量。

   需要在进程/线程之间平均分配任务从而满足条件 ①，这称为 负载均衡。

   ② 这个分配可以使得需要的 通信量是最小 的。

   并行化：将串行程序或算法转换为并行程序的过程。

2. 安排 进程/线程 之间的同步

3. 安排 进程/线程 之间的通信

2.4.3 共享内存

共享内存系统中，变量是可以 共享的(shared) 或 私有的(private)。

• 共享变量可以被任何线程读、写，而私有变量只能被单个线程访问。

• 线程间的通信是通过共享变量实现的，所以通信是隐式的。

• 动态线程和静态线程

  • 许多情况下，共享内存程序使用的是 动态线程。

    • 在这种范式中，有一个 主线程，并在任何时刻都有 一组工作线程（可能为空）。主线程通常 等待工作请求（例如，通过网络）。
    • 当一个请求到达时，它 派生出一个工作线程 来执行该请求。
    • 当工作线程 完成任务，就会 终止执行再合并到主线程中。

    这种模式充分利用了系统的资源，因为线程需要的资源只再线程实际运行时使用。

  • 另一种程序运行模式是 静态线程 范式。

    • 在这种范式中，主线程再完成必需的设置后，派生出所有的线程，在工作结束前所有的线程都在运行。
    • 当所有的线程都合并到主线程后，主线程需要做一些清理工作（如释放内存），然后也终止。

    这种范式在资源利用方面不是很高效。但是，线程的派生和合并操作是耗时的，所以如果所需的资源是可用的，静态线程模式可能获得更高的性能。它也更接近于分布式内存编程中最广泛使用的模式。

• 非确定性 (Nondeterminism)

  • 在任何一个 MIMD 系统中，如果处理器 异步执行，那么很可能会引发 非确定性。

  • 给定的输入能产生不同的输出，这种计算称为 非确定性。

  • 如果多个线程独立执行任务，每次运行时它们完成语句的速度各不相同，那么程序的结果也不同。

  示例1：线程0的 my_x 的值是 7，线程1 的 my_x 的值是 19

  C

  折叠 复制

  printf("Thread %d > my_val = %d\n", my_rank, my_x);

  两种输出可能：

  PLAINTEXT

  折叠 复制

  Thread 0 > my_val = 7
  Thread 1 > my_val = 19

  也可能是：

  PLAINTEXT

  折叠 复制

  Thread 1 > my_val = 19
  Thread 0 > my_val = 7

  • 线程独立执行并且独立地与操作系统进行交互，执行一个线程完成一段语句所花的时间在不同次的执行也是不同，所以语句执行的顺序是不能预测的。

  • 在共享内存程序中，非确定性容易导致程序错误。

  示例2：假设一个线程计算一个 int 型整数，这个 int 型整数存储在私有变量 my_val 中。假设想要将 my_val 的值加到共享内存的 x 位置中，x 初始化为 0.两个线程都要执行下面的代码：

  C

  折叠 复制

  my_val = Compute_val(my_rank);
  x += my_val;

  一次加法操作通常需要将两个数加载到寄存器、相加、最后存储结果。为了使过程简单化，假设值加载时，是从主存中直接加载到寄存器中；在存储时，直接将值从寄存器存储到主存中。

  下面是一个可能的时间顺序：

  

  显然这不是我们想要的，也很容易想象到，有其他事件顺序会产生一个 x 的不正确值。这里的非确定性是两个线程尝试同时更新内存区域 x 而造成的。

  • 当线程或者进程尝试 同时访问一个资源 时，这种访问会引发错误，我们经常说程序有 竞争条件。

    在上面的例子中，线程竞争执行 x += my_val 。这种情况下，除非一个线程在另一个线程开始前，计算完成 x += my_val ，结果才是正确的。

  • 一次只能被一个线程执行 的代码 称为 临界区。

  • 保证互斥执行的最常用机制是 互斥锁 (mutual exclusion clock) 或 互斥量 (mutex)，或者 锁 (lock)。

    互斥量 是由硬件支持的一个特殊类型的对象。基本思想是：每个临界区由一个锁来保护。在一个线程能够执行临界区的代码前，它必须通过调用一个互斥量函数来获取互斥量，在执行完临界区代码时，通过调用解锁函数来释放互斥量。

  • 监视器：在更高层次提供互斥执行。监视器是一个对象，这个对象的方法，一次只能被一个线程执行。

  • 当一个线程 “拥有” 锁时，即从调用加锁函数返回但还没有调用解锁函数时，其他线程尝试执行临界区中的代码必须在调用加锁函数时等待。

  因此，为了保证代码正常运行，我们必须修改代码，使它看起来像：

  C

  折叠 复制

  my_val = Compute_val(my_rank);
  Lock(&add_my_val_lock);
  x += my_val;
  Unlock(&add_my_val_lock);

  使用互斥量加强了临界区的串行性。

  • 替代互斥量的方案：忙等待，一个线程进入一个循环，这个循环的目的只是测试一个条件。
  C

  折叠 复制

  my_val = Compute_val(my_rank);
  if(my_rank == 1)
      while(!ok_for_1);  /* 忙等待 */
  x += my_val;
  if(my_rank == 0)
      ok_for_1 = true;

• 线程安全性

  • 谨慎使用 有 静态局部变量 的函数。函数中声明的静态变量，在函数调用时不会被销毁。因此，静态变量能够被调用函数的 线程共享，这会引起无法预测的后果。

  • C 语言 string 库函数的 strtok ，将一个输入字符串分成多个子字符串。当它第一次调用时，传递一个待分割的字符串，而在随后的调用中，它返回分割好的连续的字串。在第一次调用时，通过一个 静态的 char * 变量来指示传递给它的字符串。

    现在假设有两个线程要将字符串分割成子串。首先，线程0对 strtok 进行了首次调用，然后在 线程 0 完成分割操作完成前，线程1就对 strtok 的首次调用，这样会导致线程 0 的字符串丢失或者覆盖，在随后的调用时它会得到线程 1 的子串。

  • 类似 strtok 这样的函数不是 线程安全 的。如果它被多线程程序使用，那么会产生错误或者未知结果。当一段代码不是线程安全的，通常是因为不同的线程在访问共享的数据。

2.4.4 分布式内存

• 在分布式内存程序中，各个核能够直接访问自己的私有内存。

• 分布式内存程序通常执行多个进程。

• 消息传递

  消息传递的 API （至少）要提供一个发送和一个接收函数。进程之间通过它们的序号相互识别。序号的范围从 0~p-1，其中 p 表示进程的个数。

  例如，进程1 可以使用下面的伪代码向进程0发送消息：

  PLAINTEXT

  折叠 复制

  char message[100];
  ...;
  my_rank = Get_rank();
  if(my_rank == 1){
      sprintf(message, "Greetings from process 1");
      Send(message, MSG_CHAR, 100, 0);
  }else if(my_rank == 0){
      Receive(message, MSG_CHAR, 100, 1);
      printf("Process 0 > Received: %s\n", message);
  }

  • 程序段是 SPMD，两个进程使用相同的可执行代码，但执行不同的操作。
  • 不同进程中，变量 message 指向不同的内存块。
  • 调用 Send 函数最简单的行为是 阻塞。
  • 广播 通信：单个进程发送相同的数据给所有的进程。
  • 归约 (reduction)，在归约函数中，将各个进程计算的结果汇总成一个结果。

• 单向通信

  • 消息传递中，一个进程必须调用一个发送函数，并且发送函数必须与另一个进程调用的接收函数相匹配。任何通信都需要两个进程的显式参与。
  • 单向通信 或 远程内存访问中，单个处理器调用一个函数。在这个函数中，或者用来自另一个进程的值来更新局部内存，或者使用来自调用进程的值更新远端内存。
  • 简化通信，消除进程间同步的代价，降低通信的开销。

2.5 输入和输出

当并行程序需要输入/输出时，我们会做一些假设并遵循一些规则：

• 在分布式内存程序中，只有进程 0 能够访问 stdin。

  在共享内存程序中，只有主线程或线程0能够访问 stdin

• 在分布式内存和共享内存系统中，所有线程/进程都能够访问 stdout 和 stderr 。

• 因为输出到 stdout 的非确定性顺序，大多数情况下，只有一个进程/线程会将结果输出到 stdout 。但输出调试程序的结果是一个例外。这种情况下，允许多个进程/线程写 stdout 。

• 只有一个进程/线程会尝试访问一个除 stdin\stdout 或者 stderr 外的文件。所以，每个进程/线程能够打开自己的私有文件进行读、写，但是没有两个进程/线程能够打开相同的文件。

• 调试程序输出在生成输出结果时，应该包括进程/线程的序号或者进程标识符。

2.6 性能

2.6.1 加速比和效率 ⭐

当在 p 核系统上运行程序，每个核运行一个进程或线程，并行程序的运行速度就是串行程序速度的 p 倍。如果我们称 串行运行时间为 串行𝑇串行 ，并行运行时间为 并行𝑇并行，那么最佳的预期是 并行串行𝑇并行=𝑇串行/𝑝 ，此时，我们称并行程序有 线性加速比。

程序的 加速比：

串行并行𝑆=𝑇串行𝑇并行

线性加速比为 S = p，这非常难达到。

效率：

串行并行串行并行𝐸=𝑆𝑝=𝑇串行𝑇并行𝑝=𝑇串行𝑝⋅𝑇并行

如果用 开销𝑇开销 表示并行开销，那么

并行串行开销𝑇并行=𝑇串行/𝑝+𝑇开销

2.6.2 阿姆达尔定律

• 阿姆达尔定律 描述了这样的一个事实：大致上，除非一个串行程序的执行几乎全部都并行化，否则，无论有多少可以利用的核，通过并行化产生的加速比都是受限的。

示例：一个串行程序中，可以并行化其中的 90%。假设如果使用 𝑝 个核，则程序可并行化部分的加速比就是 𝑝 。若该程序串行版本运行时间 串行秒𝑇串行=20秒 。

并行化后，其中的可并行部分的运行时间就是：串行秒90%×𝑇串行/𝑝=18/𝑝秒

不可并行化部分的运行时间为 串行秒10%×𝑇串行=2秒

加速比为：𝑆=2018/𝑝+2⩽10

2.6.3 可扩展性

• 假设运行一个拥有固定进程或线程数目的并行程序，并且它的输入规模也是固定的，那么我们可以得到一个效率值 E。
• 增加该程序所用的进程/线程数，如果在输入规模也以相应增长率增加的情况下，该程序的效率值一直都是 E，那么我们就称该程序是 可扩展的。

示例：程序串行版本的运行时间是 串行𝑇串行=𝑛，其中，串行𝑇串行 的单位是微秒，𝑛 可以看作是问题规模。假设程序并行版本的运行时间是 串行𝑇串行=𝑛/𝑝+1 。那么，效率值就为：

𝐸=𝑛𝑝(𝑛/𝑝+1)=𝑛𝑛+𝑝

如果程序是可扩展的，以 𝑘 为倍率提高进程/线程的数目，希望在效率值 𝐸 不变的情况下，找到问题规模的增加比例 𝑥。增加进程/线程的个数，使之变为 𝑘𝑝 ，问题的规模相应增加到 𝑥𝑛 ，通过求解方程来求𝑥 的值：

𝐸=𝑛𝑛+𝑝=𝑥𝑛𝑥𝑛+𝑘𝑝

如果 𝑥=𝑘 则等式成立。

只有问题规模增加的倍率与进程/线程数增加的倍率相同时，效率才会是恒定的，从而程序是可扩展的。

• 如果在增加进程/线程的个数时，可以维持固定的效率，却不增加问题的规模，那么程序称为 强可扩展性的。
• 如果在增加进程/线程个数的同时，只有以相同倍率增加问题的规模才能使效率值保持不变，那么程序就称为 弱可扩展的 (weakly scalable)。

2.6.4 计时

• 标准 C 函数 clock 所提高的时间，是 CPU 时间，代表程序执行代码的总时间，不包括程序空闲状态的时间。
• 我们关心：代码从开始执行到执行结束的总耗费时间。

double start, finish;
...;
start = Get_current_time();
...;
finish = Get_current_time();
printf("The elapsed time = %e seconds\n", finish-start);

• MPI：MPI_Wtime
• OpenMP：omp_get_wtime

并行程序中，计时的代码可能被多个进程/线程执行，这将导致最后输出 𝑝 个运行时间。我们通常感兴趣的是单个时间：从第一个进程/线程开始执行，到最后一个进程/线程结束之间所花费的时间。

shared double global_elapsed;
private double my_start, my_finish, my_elapsed;
...;
/* barrier 函数：同步所有进程/线程 */
Barrier();
my_start = Get_current_time();

/* 需要计时的执行代码 */
...;

my_finish = Get_current_time();
my_elapsed = my_finish - my_start;

/* 在所有进程/线程中找最大 */
global_elapsed = Global_max(my_elapsed);
if(my_rank == 0)
    printf("The elapsed time = %e seconds\n", global_elapsed);

2.7 并行程序设计

串行程序并行化：（Foster 方法）

1. 划分 (partitioning)：将要执行的指令核数据按照计算部分拆分成多个小任务。这一步的关键在于识别出可以并行执行的任务。
2. 通信 (communication)：确定上一步所识别出来的任务之间需要执行那些通信。
3. 凝聚或聚合 (agglomeration or aggregation)：将第一步所确定的任务与通信结合成更大的任务。
4. 分配 (mapping)：将上一步聚合好的任务分配到进程/线程中。这一步还要使通信量最小化，使各个进程/线程所得到的工作量大致均衡。

题目：

2.10 假设一个程序需要运行 1012 条指令来解决一个特定问题，一个单处理器系统可以在 106 秒（大约11.6天）内完成。所以，一个单处理器系统平均每秒运行 106 条指令。现在假设程序已经实现并行化，可以在分布式内存系统上运行。 该并行程序使用 p 个处理器，每个处理器执行 1012/𝑝 条指令并必须发送 109(𝑝−1) 条消息。执行该程序时，不会有额外的开销，即每个处理器执行完所有的指令并发送完所有的消息之后，程序就完成了，而不会有由诸如等待消息等事件所产生的延迟。那么， a. 假设发送一条消息需要耗费 10−9 秒。如果程序使用1000个处理器，每个处理器的速度和单个处理器运行串行程序的速度一样，那么该程序的运行需要多少时间？ b. 假设发送一条消息需要耗费 10−3 秒。如果程序使用1000个处理器，那么该程序的运行需要多少时间？

a. 秒1012/10001012/1000∗109∗999∗10−9+106/(1000)=1999秒

b. 秒1012/10001012/1000∗109∗999∗10−3+106/(1000)=999001000秒

3 用 MPI 进行分布式内存编程

消息传递接口(Message-Passing Interface, MPI)

3.1 预备知识

经典程序：

#include <stdio.h>
int main(void){
    printf("hello, world\n");
    
    return 0;
}

使用 MPI 编写一个类似的 “hello, world” 程序，不让每个进程都简单地打印一条消息，相反，指派其中一个进程负责输出，其他进程向它发送要打印的消息。

指令 0号进程为输出进程，其余进程向它发送消息。

3.1.1 编译与执行

编译：

mpicc -g -Wall -o mpi_hello mpi_hello.c

打印来自进程问候句的 MPI 程序：

#include <stdio.h>
#include <string.h>
#include <mpi.h>

const int MAX_STRING = 100;

int main(void){
    char greeting[MAX_STRING];
    int  comm_sz;  /* 进程数量 */
    int  my_rank;  /* 我的进程编号 */
    
    MPI_Init(NULL, NULL);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    
    if(my_rank != 0){
        sprintf(greeting, "Greeting from process %d of %d!",
               my_rank, comm_sz);
        MPI_Send(greeting, strlen(greeting) + 1, MPI_CHAR, 0, 0,
                MPI_COMM_WORLD);
    } else {
        printf("Greetings from process %d of %d!\n", my_rank,
               comm_sz);
        for(int q = 1; q < comm_sz; q++){
            MPI_Recv(greeting, MAX_STRING, MPI_CHAR, q, 0, 
                     MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            printf("%s\n", greeting);
        }
    }
    MPI_Finalize();
    return 0;
} /* main */

启动程序：

mpiexec -n <number of processes> ./mpi_hello

mpiexec -n 4 ./mpi_hello

3.1.2 MPI 程序

一个 C 语言程序，包含了 C 语言的标准头文件。

包含 mpi.h 头文件。

3.1.3 MPI_Init 和 MPI_Finalize

MPI_Init ：告知 MPI 系统进行必要的初始化设置。

❗在调用 MPI_Init 之前不应该调用其他 MPI 函数。

int MPI_Init(
    int*    argc_p,  /* in/out */
    char*** argv_p   /* in/out */
)

• 参数 argc_p 和 argv_p 是指向参数 argc 和 argv 的指针。程序不需要使用时设置 NULL。
• 返回一个 int 型错误码。

MPI_Finalize ：告知 MPI 系统 MPI 已经执行完毕，为 MPI 而分配的任何资源都可以释放了。

int MPI_Finalize(void);

❗在调用 MPI_Finalize 之后不应该调用其他 MPI 函数。

MPI 程序的基本框架：

...
#include <mpi.h>
...

int main(int argc, char* argv[]){
    ...;
    /* 此前没有任何 MPI 函数调用 */
    MPI_Init(&argc, &argv);
    ...;
    MPI_Finalize();
    /* 此后没有任何 MPI 函数调用 */
    return 0;
}

3.1.4 通信子、MPI_Comm_size 和 MPI_Comm_rank

通信子 (communicator) ：一组可以互相发送消息的进程集合。

MPI_Init 的其中一个目的：在用户启动程序时，定义由用户启动的所有进程所组成的通信子，这个通信子称为 MPI_COMM_WORLD 。

int MPI_Comm_size(
    MPI_Comm comm       /* in  - 通信子 */,
    int*     comm_sz_p  /* out - 进程的数量*/
);

int MPI_Comm_rank(
    MPI_Comm comm       /* in  - 通信子 */,
    int*     my_rank_p  /* out - 进程号 */
);

3.1.5 SPMD 程序

单程序多数据流 (Single Program, Multiple Data, SPMD)：编写一个单个程序，让不同进程产生不同的动作。实现的方式是，简单地让进程按照它们的进程号来匹配程序分支。

3.1.6 通信

0 号进程负责打印输出消息，使用一个 for 循环接收并打印由 1、2、···、comm_sz-1 号进程发送来的消息。

3.1.7 MPI_Send

int MPI_Send(
    void*         msg_buf_p     /* in - 包含消息内容的内存块指针*/,
    int           msg_size      /* in - 数据大小*/,
    MPI_Datatype  msg_type      /* in - 数据类型*/,
    int           dest          /* in - 接收消息的进程号*/,
    int           tag           /* in - 区分看上去完全一样的消息*/,
    MPI_Comm      communicator  /* in - 通信子*/
);

• 消息长度：

  字符串大小为 strlen(str) + 1

• 数据类型：

  MPI 数据类型	C 语言数据类型
  MPI_CHAR	char
  MPI_SHORT	short
  MPI_INT	int
  MPI_LONG	long
  MPI_LONG_LONG	long long
  MPI_UNSIGNED_CHAR	unsigned char
  MPI_UNSIGNED_SHORT	unsigned short
  MPI_UNSIGNED	unsigned int
  MPI_UNSIGNED_LONG	unsigned long
  MPI_UNSIGNED_LONG_LONG	unsigned long long
  MPI_FLOAT	float
  MPI_DOUBLE	double
  MPI_LONG_DOUBLE	long double
  MPI_BYTE	unsigned char
  MPI_PACKED	(用于打包数据)

• tag 可以用于一个进程向另一个进程发送数据，其中一部分数据用于输出，可以设置 tag = 0，另一部分数据用于计算，可以设置 tag = 1，以让接收进程区分数据的使用。

• 通信子重要作用功能：指定通信范围。

3.1.8 MPI_Recv

int MPI_Recv(
    void*        msg_buf_p    /* out - 内存块指针*/, 
    int          buf_size     /* in  - 要存储对象的数量*/,
    MPI_Datatype buf_type     /* in  - 对象的类型*/,
    int          sourse       /* in  - 接收消息的来源进程号*/,
    int          tag          /* in  - 要与发送消息的 tag 参数匹配*/,
    MPI_Comm     communicator /* in  - 与发送消息的通信子参数匹配*/,
    MPI_Status*  status_p     /* out*/
);

status_p 一般不使用这个参数，赋予 MPI_STATUS_IGORE 就行。

3.1.9 消息匹配

• 假定 q 号进程调用了 MPI_Send 函数：

  C

  折叠 复制

  MPI_Send(send_buf_p, send_buf_sz, send_type, dest, send_tag,
          send_comm);

• 假定 r 号进程调用了 MPI_Recv 函数：

  C

  折叠 复制

  MPI_Recv(recv_buf_p, recv_buf_sz, recv_type, serc, recv_tag, 
          recv_comm, &status);

• 则 q 号进程调用 MPI_Send 函数所发送的消息可以被 r 号进程调用 MPI_Recv 函数接收，如果：

  • recv_comm = recv_comm

  • recv_tag = send_tag

  • dest = r

  • src = q

  • 前三对参数：send_buf_p/recv_buf_q、send_buf_sz/recv_buf_sz 和 send_type/recv_type 必须指定兼容的缓冲区。

    大多数情况，满足下面条件即可：如果 recv_type = send_type ，同时 recv_buf_sz ≥ send_buf_sz ，那么由 q 号进程发送的消息就可以被 r 号进程成功地被接收。

MPI 提供了一个特殊的常量 MPI_ANY_SOURCE，可以传递给 MPI_Recv 的参数 source ，这样，如果 0 号进程执行下列代码，它可以按照进程完成工作的顺序来接收结果：

for(i = 1; i<comm_sz; i++){
    MPI_Recv(result, result_sz, result_type, MPI_ANY_SCOURCE, 
            result_tag, comm, MPI_STATUS_IGNORE);
    Process_result(result);
}

类似地，一个进程可能接收多条来自另一个进程的有着不同标签的消息，并且接收进程并不知道消息发送的顺序。在这种情况下，MPI 提供了特殊的常量 MPI_ANY_TAG ，将它传给 MPI_Recv 的参数 tag 。

3.1.10 status_p 参数

接收者可以在不知道以下信息的情况下接收消息：

1. 消息的数据量
2. 消息的发送者
3. 消息的标签

接收者如何找出这些值的？MPI_Recv 的最后一个参数的类型为 MPI_Status * ，MPI_Status 是一个由至少三个成员的结构：MPI_SOURCE、MPI_TAG 和 MPI_ERROR 。假定程序含有如下定义：

MPI_Status status;

那么将 &status 作为最后一个参数传递给 MPI_Recv 函数并调用它后，可以通过检查以下两个成员来确定发送者和标签：

status.MPI_SOURCE;
status.MPI_TAG;

接受到的数据量不是存储在应用程序可以直接访问的域中，但用户可以调用 MPI_Get_count 函数找回这个值。例如对 MPI_Recv 的调用中，接收缓冲区的类型为 recv_type，再次传递 &status 参数，则以下调用

MPI_Get_count(&status, recv_type, &count);

会返回 count 参数接收到的元素数量。

MPI_Get_count 的语法结构：

int MPI_Get_count(
    MPI_Status*  status_p /* in */,
    MPI_Datatype type     /* in */,
    int*         count_p  /* out */
);

3.1.11 MPI_Send 和 MPI_Recv 的语义

将消息从一个进程发送到另一个进程：

• 发送进程组装消息。

• 发送进程可以 缓存 消息也可以 阻塞。

  • 如果它缓冲消息，则 MPI 系统将会把消息放置在它自己的内部存储器里，并返回 MPI_Send 的调用。
  • 如果发生阻塞，那么它将一直等待，直到可以发送消息，并不立即返回 MPI_Send 的调用。

  如果使用 MPI_Send，当函数返回时，实际上并不知道消息是否已经发送出去，我们只知道消息所用的存储器，即发送缓冲区，可以被程序再次使用。

• MPI_Send 的精确行为由 MPI 实现所决定，典型的实现方法有一个 “截止”大小（“cutoff” message size）。如果一个消息的大小小于“截止”大小，它将被缓冲；如果大于 “截止”大小，那么 MPI_Send 将会被阻塞。

• MPI_Recv 函数总是阻塞的，知道接收到一条匹配的消息。因此，当 MPI_Recv 函数返回时，就知道一条消息已经存储在接收缓冲区中（除非产生了错误）。

• 非阻塞的 MPI_lsend 和 MPI_lrecv 会立即返回。

3.1.12 潜在的陷阱

• ❗MPI_Recv 的语义会导致 MPI 编程中的一个潜在陷阱：如果一个进程试图接收消息，但没有相匹配的消息，那么该进程将会被永远阻塞在那里，即 进程悬挂。
• ❗如果调用 MPI_Send 发生了阻塞，并且没有相匹配的接收，那么发送进程就悬挂起来，如果调用 MPI_Send 被缓冲，但没有相匹配的接收，那么消息将被丢失。

3.2 用 MPI 来实现梯形积分法

梯形面积和 = ℎ[𝑓(𝑥0)/2+𝑓(𝑥1)+⋯+𝑓(𝑥𝑛−1)+𝑓(𝑥𝑛)/2]

串行程序的伪代码：

h = (b-a)/n;
approx = (f(a) + f(b))/2.0;
for(i=1; i<=n-1; i++){
    x_i = a + i*h;
    approx += f(x_i);
}
approx = h*approx;

并行化梯形积分法

四个基本步骤设计一个并行程序：

1. 将问题的解决方案划分成多个任务
2. 在任务间识别出需要的通信
3. 将任务聚合成复合任务
4. 在核上分配复合任务

识别两种任务：一是获取单个矩形区域的面积，另一种是计算这些区域的面积和。然后利用通信信道将第一个任务和第二个任务相连接。

将梯形区域面积的计算聚合成组。实现这一目标：将区间 [a, b] 分成 comm_sz 个子区间。如果 comm_sz 可以整除 n ，我们可以简单地在 n/comm_sz 个梯形和所有 comm_sz 个子空间上应用梯形积分法。最后，我们可以利用进程种的某一个，如 0 号进程，我们将这些梯形面积的估计值累加起来完成整个计算过程。

Get a, b, n;
h = (b-a)/n;
local_n = n/comm_sz;
local_a = a + my_rank * local_n * h;
local_b = local_a + local_n * h;
local_integral = Trap(loca_a, local_b, local_n, h);

if(my_rank != 0)
    发送 local_integral 给 进程0;
else{ /* my_rank == 0 */
    total_integral = local_integral;
    for(proc = 1; proc < comm_sz; proc++){
        从进程 proc 种接收 local_integral;
        total_integral += local_integral;
    }
}
if(my_rank == 0)
    print result;

3.3 I/O 处理

3.3.1 输出

如果多个进程试图写标准输出 stdout，那么这些进程的输出顺序是无法预测的，甚至会发送一个进程的输出被另一个进程打断的情况。

这是因为，MPI 进程都在相互“竞争”，以取得共享输出设备、标准输出 stdout 的访问。这种竞争导致不确定性。

梯形积分法 MPI 程序的第一个版本：

int main(void){
    int my_rank, comm_sz, n = 1024, local_n;
    double a = 0.0, b = 3.0, h, local_a, local_b;
    double local_int, total_int;
    int source;
    
    MPI_Init(NULL, NULL);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    
    h = (b-a)/n;
    local_n = n/comm_sz;
    
    local_a = a + my_rank * h * local_n;
    local_b = local_a + h * local_n;
    local_int = Trap(local_a, local_b, local_n, h);
    
    if(my_rank != 0){
        MPI_Send(&local_int, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
    }else{
        total_int = local_int;
        for(int source = 0; source<comm_sz; source++){
            MPI_Recv(&local_int, 1, MPI_DOUBLE, source, 0, 
                     MPI_COMM_WORLD, MPI_STATUS_IGNORE);
            total_int += local_int;
        }
    }
    
    if(my_rank == 0){
        printf("With n = %d trapezoids, our estimate\n", n);
        printf("of the integral from %f to %f = %.15e\n", 
               a, b, total_int);
    }
    
    MPI_Finalize();
    return 0;
} /* main */

梯形积分法里的 Trap 函数：

double Trap(
    double left_endpt, 
    double right_endpt,
    int trap_count,
    double base_len
){
    double estimate, x;
    int i;
    
    estimate = (f(left_endpt) + f(right_endpt))/2.0;
    for(int i = 1; i <= trap_count, i++){
        x = left_endpt + i*base_len;
        estimate += f(x);
    }
    estimate = estimate * base_len;
    
    return estimate;
}

3.3.2 输入

大部分的 MPI 实现只允许 MPI_COMM_WORLD 中的 0 号进程访问标准输入 stdin 。

根据进程号选取转移分支。0 号进程负责读取数据，并将数据发送给其他进程。

梯形积分法的并行程序中的 Get_input 函数：

void Get_input(
    int     my_rank   /* in */,
    int     comm_sz   /* in */,
    double* a_p       /* out */,
    double* b_p       /* out */,
    int*    n_p       /* out */
){
    int dest;
    
    if(my_rank == 0){
        printf("Enter a, b and n\n");
        scanf("%lf %lf %d", a_p, b_p, n_p);
        for(dest = 1; dest < comm_sz; dest++){
            MPI_Send(a_p, 1, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD);
            MPI_Send(b_p, 1, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD);
            MPI_Send(n_p, 1, MPI_INT, dest, 0, MPI_COMM_WORLD);
        }
    }else{
        MPI_Recv(a_p, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD,
                 MPI_STATUS_IGNORE);
        MPI_Recv(b_p, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD,
                 MPI_STATUS_IGNORE);
        MPI_Recv(n_p, 1, MPI_INT, 0, 0, MPI_COMM_WORLD,
                 MPI_STATUS_IGNORE);
    }
} /* Get_input */

❗我们必须在初始化 my_rank 和 comm_sz 后，才能调用该函数

...;
MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

Get_input(my_rank, comm_sz, &a, &b, &c);

3.4 集合通信

梯形积分法程序中，“全局求和” 部分，0 号进程要做全局求和的所有工作。希望有一个更公平的分配方法。

3.4.1 树形结构通信

可以使用 二叉树结构，一开始 1 号进程、3 号进程、5 号进程、7号进程将它们的值发送给 0 号进程、2 号进程、4 号进程、6 号进程。然后 0 号进程、2 号进程、4 号进程、6 号进程将接收到的值加到它们自己原有的值上，整个过程重复两次：

1. ① 2 号进程和6号进程将它们的新值分别发送给 0号进程、4 号进程。

   ② 0 号进程和4号进程将接收到的值加到它们自己新值上

2. ① 4 号进程将它们的新值分别发送给 0号进程

   ② 0 号进程将接收到的值加到它们自己新值上

新方案：很多工作是由不同进程并发完成。例如，如何 comm_sz = 1024，则原方案需要 0 号进程 执行 1023 次接收和加法操作，新方案中，0号进程只需要 10 次接收和加法操作。

3.4.2 MPI_Reduce

“全局求和函数” 需要通信。全局求和函数可能涉及两个及以上的进程。在 梯形积分法程序中，它涉及 MPI_COMM_WORLD 中所有的进程。

• 在 MPI 里，涉及 通信子中所有进程 的通信函数称为 集合通信 (collective communication)。
• MPI_Send 和 MPI_Recv 通常称为 点对点通信 (point-to-point communication)

全局求和函数只是集合通信函数类别中的一个特殊例子而已。例如，我们可能并不是想计算分布在各个进程上的 comm_sz 值的总和，而是想知道最大值和最小值，或者总的乘积，或者其他许多可能情况中的任何一个产生的结果。

MPI 对全局求和函数进行概括，使这些可能性中的任意一种用单个函数实现：

int MPI_Reduce(
    void*        input_data_p   /* in */,
    void*        output_data_p  /* out */,
    int          count          /* in */,
    MPI_Datatype datatype       /* in */,
    MPI_Op       operator       /* in */,
    int          dest_process   /* in */,
    MPI_Comm     comm
);

使用 MPI_Reduce 求和：

MPI_Reduce(&local_int, &total_int, 1, MPI_DOUBLE, MPI_SUM, 0, 
           MPI_COMM_WORLD);

如果 count 参数大于 1，那么 MPI_Reduce 函数可以应用到数组上，而不仅仅是应用在简单的标量上。下面的代码可以用于 一组 N 维向量的加法，每个进程上有一个向量：

double local_x[N], sum[N];
...;
MPI_Reduce(local_x, sum, N, MPI_DOUBLE, MPI_SUM, 0, 
           MPI_COMMM_WORLD);

MPI 中预定义的归约操作符：

3.4.3 集合通信与点对点通信

1. 在 通信子 中的 所有进程 必须调用相同的集合通信函数。

2. 每个进程传递给 MPI 集合通信函数的参数必须是相容的。

   ❗如果一个进程将 0 作为 dest_process 的值传递给函数，另一个传递的是 1 ，那么对 MPI_Reduce 调用就是错误的。

3. 参数 output_data_p 只用在 dest_process 上，然而，所有进程仍需传递一个与 output_data_p 相对应的实际参数，即使它的值只是 NULL。

4. 点对点通信函数 是通过标签和通信子来匹配的。

   ❗集合通信函数 不使用标签，只通过 通信子 和 调用的顺序 来进行匹配。

   示例：

   

   假设上面 MPI_Reduce 运算符均为 MPI_SUM，目标进程都是 0 号进程。两次调用 MPI_Reduce 后，0 号进程的 b = 1+2+1 = 4，d = 2 + 1 + 2 = 5。

   ❗不能使用 同一个缓冲区 同时作为 输入和输出调用 MPI_Reduce ，下面这种调用方式是非法的，它的结果不可预测，它可能产生一个错误结果，也可能导致程序崩溃，但也可能产生正确的结果：

   C

   折叠 复制

   /*❌非法代码*/MPI_Reduce(&x, &x, 1, MPI_DOUBLE, MPI_SUM, 0, com);

3.4.4 MPI_Allreduce

用一棵树来计算全局总和，我们可以“颠倒” (reverse) 整棵树来发布全局总和。

还有另一种替代方法，可以让进程之间相互交换部分结果，而不是单向的通信。这种通信模式称为 蝶形：

MPI 提供 MPI_Reduce 的变种—— MPI_Allreduce，可以令通信子中的 所有进程都存储结果：

int MPI_Allreduce(
    void*        input_data_p   /* in */,
    void*        output_data_p  /* out */,
    int          count          /* in */,
    MPI_Datatype datatype       /* in */,
    MPI_Op       operator       /* in */,
    MPI_Comm     comm           /* in */
)

参数表与 MPI_Reduce 相同，除了没有 dest_process 这个参数，因为所有进程都能得到结果。

3.4.5 广播

简单地在树形全局求和里“颠倒”通信，我们可以得到下面树形结构通信图，并且能够将这个结构用于输入数据的发布。

广播：属于一个进程的数据被发送到通信子中的所有进程 的集合通信。

int MPI_Bcast(
    void*        data_p      /* in/out❗ */,
    int          count       /* in */,
    MPI_Datatype datatype    /* in */,
    int          source_proc /* in - 源进程号 */,
    MPI_Comm     comm        /* in */ 
);

进程号为 scourse_proc 的进程将 data_p 所引用的内存内容发送给了通信子 comm 的所有进程。

使用 MPI_Bcast 的 Get_input 函数：

void Get_input(
    int     my_rank   /* in */,
    int     comm_sz   /* in */,
    double* a_p       /* out */,
    double* b_p       /* out */,
    int*    n_p       /* out */
){
    int dest;
    
    if(my_rank == 0){
        printf("Enter a, b and n\n");
        scanf("%lf %lf %d", a_p, b_p, n_p);
    }
    MPI_Bcast(a_p, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(b_p, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(n_p, 1, MPI_INT, 0, MPI_COMM_WORLD);
} /* Get_input */

3.4.6 数据分发

编写一个程序，用于计算向量和：

𝑥+𝑦=(𝑥0,𝑥1,⋯,𝑥𝑛−1)+(𝑦0,𝑦1,⋯,𝑦𝑛−1)=(𝑥0+𝑦0,𝑥1+𝑦1,⋯,𝑥𝑛−1+𝑦𝑛−1)=(𝑧0,𝑧1,⋯,𝑧𝑛−1)=𝑧

如果用 double 类型的数组来表示向量，可以用串行加法来计算向量求和：

void Vector_sum(double x[], double y[], double z[], int n){
    int i;
    
    for(i = 0; i<n; i++){
        z[i] = x[i] + y[i];
    }
}  /* Vector_sum */

如何用 MPI 来实现这个程序？

• 指定各个任务求和的对应分量。各个任务间没有通信，向量并行加法问题归结为聚合任务以及将它们分配到核上。
• 如果分量的个数为 n，并且我们有 comm_sz 个核或者进程，那么可以简单地将连续 local_n 个向量分量所构成的块，分配到每个进程中，这种做法称为 向量 的 块划分。
• 循环划分：轮转方式分配向量分量。
• 块-循环划分：用一个循环分发 向量分量所构成的块，而不是分发单个向量分量。

无论哪种划分方法，每个进程都有一个 local_n 个向量分量，用一个只有 local_n 个元素的数组存储这些分量:

void Parallel_vector_sum(
    double local_x[] /* in */,
    double local_y[] /* in */,
    double local_z[] /* out */,
    int    local_n   /* in */
){
    int local_i;
    
    for(local_i = 0; local_i < local_n; local_i++)
        local_z[local_i] = local_x[local_i] + local_y[local_i];
} /* Parallel_vector_sum */

3.4.7 散射

测试向量加法函数。首先读取向量的维度，然后读取向量 x 和 向量 y 。0 号进程读取输入向量，然后将它们广播给其他进程，但这种方法很浪费。每个进程只要读取它所需要操作的部分分量。

因此，可以 0 号进程读入整个向量，但只将分量发送给需要分量的其他进程。MPI 提供这样的函数：

int MPI_Scatter(
    void*        send_buf_p  /* in */,
    int          send_count  /* in - 发送到每个进程的数据量*/,
    MPI_Datatype send_type   /* in */,
    void*        recv_buf_p  /* out */,
    int          recv_count  /* in - 每个进程接收的数据量*/,
    MPI_Datatype recv_type   /* in */,
    int          src_proc    /* in - 源进程号 */,
    MPI_Comm     comm        /* in */
);

如果 通信子 comm 包含 comm_sz 个进程，那么 MPI_Scatter 函数会将 send_buf_p 所引用的数据分成 comm_sz 份，第一份给 0 号进程，第二份给 1 号进程，第三份给 2 号进程，以此类推。

例如，假如我们使用块划分法，并且 0 号进程已经将一个有 n 个分量的向量整个读入 send_buf_p 中，则 0 号进程将得到第一组 local_n = n/comm_sz 个分量，1 号进程将得到下一个 local_n 个分量，以此类推。每个进程应该将它本地的向量作为 recv_buf_p 参数的值，将 local_n 作为 recv_count 参数的值，send_type 和 recv_type 参数的值都应该是 MPI_DOUBLE，src_proc 参数的值应该是 0。

send_count 参数的值也应该是 local_n ，因为 send_count 表示 发送到每个进程的数据量。

注意：MPI_Scatter 函数将 send_count 个对象所组成的第一个块发送给 0 号进程，将下一个由 send_count 个对象所组成的块发送给了 1 号进程，以此类推。

所以，这种读取和分发输入向量的方法 只适用于 块划分法。

void Read_vector(
    double   local_a[]  /* out */,
    int      local_n    /* in */,
    int      n          /* in */,
    char     vec_name[] /* in */,
    int      my_rank    /* in */,
    MPI_Comm comm       /* in */
){
    double* a = NULL;
    int i;
    
    if(my_rank == 0){
        a = malloc(n * sizeof(double));
        printf("Enter the vector %s\n", vec_name);
        for(i = 0; i<n; i++)
            scanf("%lf", &a[i]);
        MPI_Scatter(a, local_n, MPI_DOUBLE, local_a, local_n,
                   MPI_DOUBLE, 0, comm);
        free(a);
    } else {
        MPI_Scatter(a, local_n, MPI_DOUBLE, local_a, local_n,
                   MPI_DOUBLE, 0, comm);
    }
} /* Read_vector */

3.4.8 聚集

我们需要编写一个可以打印分布式向量的函数。这个函数将向量的所有分量都收集到 0 号进程上，然后由 0 号进程将所有分量都打印出来，这个函数的通信由 MPI_Gather 来执行：

int MPI_Gather(
    void*        send_buf_p  /* in */,
    int          send_count  /* in */,
    MPI_Datatype send_type   /* in */,
    void*        recv_buf_p  /* out */,
    int          recv_count  /* in */,
    MPI_Datatype recv_type   /* in */,
    int          dest_proc   /* in - 目的进程号*/,
    MPI_Comm     comm
);

在 0 号进程中，由 send_buf_p 所引用的内存区的数据存储在 recv_buf_p 的第一个块中，在 1 号进程中，由 send_buf_p 所引用的内存区的数据存储在 recv_buf_p 的第二个块中，以此类推。

void Print_vector(
    double   local_b[]  /* in */,
    int      local_n    /* in */,
    int      n          /* in */,
    char     title[]    /* in */,
    int      my_rank    /* in */,
    MPI_Comm comm       /* in */
){
    double* b = NULL;
    int i;
    
    if(my_rank == 0){
        b = malloc(n * sizeof(double));
        MPI_Gather(local_b, local_n, MPI_DOUBLE, b, local_n,
                  MPI_DOUBLE, 0, comm);
        printf("%s\n", title);
        for(i = 0; i<n; i++)
            printf("%f", b[i]);
        printf("\n");
        free(b);
    } else {
        MPI_Gather(local_b, local_n, MPI_DOUBLE, b, local_n,
                  MPI_DOUBLE, 0, comm);
    }
} /* Print_vector */ 

3.4.9 全局聚集

完成矩阵和向量的相乘，如果 𝐴=(𝑎𝑖𝑗) 是一个 𝑚×𝑛 的矩阵，𝑥 是一个具有 𝑛 个分量的向量，那么 𝑦=𝐴𝑥 就是一个有 m 个分量的向量。我们可以用 A 的第 i 行与 𝑥 的点积来取 𝑦 的第 i 个分量。串行伪代码：

for(i = 0; i<m; i++){
    y[i] = 0.0;
    for(j = 0; j<n; j++)
        y[i] += A[i][j] * x[j];
}

C 语言串行函数：

void Mat_vect_mult(
    double A[]  /* in */, 
    double x[]  /* in */,
    double y[]  /* out */, 
    int m       /* in */, 
    int n       /* in */)
{
    int i, j;

    for (i = 0; i < m; i++) {
        y[i] = 0.0;
        for (j = 0; j < n; j++)
            y[i] += A[i*n+j]*x[j];
    }
} /* Mat_vect_mult */

如果 x 有一个块划分，我们如何安排使得下面循环前就能使每个进程都能访问 x 中的所有分量呢？

for (j = 0; j < n; j++)
    y[i] += A[i*n+j]*x[j];

MPI 提供了这样的一个单独的函数：

int MPI_Allgather(
    void*        send_buf_p  /* in */,
    int          send_count  /* in */,
    MPI_Datatype send_type   /* in */,
    void*        recv_buf_p  /* out */,
    int          recv_count  /* in */,
    MPI_Datatype recv_type   /* in */,
    MPI_Comm     comm        /* in */
);

void Mat_vect_mult(
    double local_A[]  /* in */,
    double local_x[]  /* in */,
    double local_y[]  /* out */,
    int local_m       /* in */,
    int n             /* in */,
    int local_n       /* in */,
    MPI_Comm comm     /* in */
) {
    double* x;
    int local_i, j;
    int local_ok = 1;
    
    x = malloc(n* sizeof(double));
    MPI_Allgather(local_x, local_n, MPI_DOUBLE, 
                 x, local_n, MPI_DOUBLE, comm);
    for(local_i = 0; local_i < m; local_i++){
        local_y[local_i] = 0.0;
        for(j = 0; j<n; j++)
            local_y[local_i] += local_A[local_i * n + j] * x[j];
    }
    
    free(x);
} /* Mat_vect_mult */

3.6 MPI 程序的性能估计

3.6.1 计时

MPI 提供的函数：

double MPI_Wtime(void);

计时：

double start, finish;
...;
start = MPI_Wtime();
...;
finish = MPI_Wtime();
printf("Proc %d > Elapsed time = %e seconds\n", 
       my_rank, finish-start);

3.6.2 结果

• 串行𝑇串行 代表串行时间，由于它取决于输入值 n，所以把它设为 串行𝑇串行(𝑛) 。同理，并行运行时间 并行𝑇并行 取决于输入值 n 和进程数目 comm_sz = p ，则设为 并行𝑇并行(𝑛,𝑝) 。并行程序会将串行程序的工作分配到各个进程上，但又会增加额外的开销，设此开销为：

  并行串行开销𝑇并行(𝑛,𝑝)=𝑇串行(𝑛)/𝑝+𝑇开销

• 在 MPI 程序中，并行计算的开销一般来自于 通信，它同时还受到问题集的规模和进程数的影响。

• 我们的矩阵-向量乘法程序并不难实现，主要的串行运算部分是一对 for 循环的嵌套：

  C

  折叠 复制

  for(i = 0; i<m; i++){
      y[i] = 0.0;
      for(j = 0; j<n; j++){
          y[i] += A[i*n+j]*x[j];
      }
  }

  如果只对浮点运算进行计数，则内层循环会执行 𝑛 次乘法和 𝑛 次加法，总共 2𝑛 次浮点运算。因为执行了 𝑚 次内层循环，则执行上面这段代码总共需要 2𝑚𝑛 次运算。所以当 𝑚=𝑛 时：串行𝑇串行(𝑛)≈𝑎𝑛2，𝑎 是常数（符号≈意味着近似等于）。

  如果串行程序要执行 𝑛×𝑛 矩阵与一个 𝑛 维向量相乘，那么并行程序中每个进程进行 𝑛/𝑝×𝑛 矩阵与 𝑛 维向量相乘。每一个局部的矩阵 - 向量相乘要执行 𝑛2/𝑝 次浮点运算，即每个进程的工作量被削减了 𝑝 倍。

  然而，并行程序在进行本地矩阵 - 向量乘法前，需要调用 MPI_Allgather 函数。在我们的例子中：

  并行串行𝑇并行(𝑛,𝑝)=𝑇串行(𝑛)/𝑝+𝑇allgather

• 𝑛 值较小，𝑝 值较大，并行𝑇并行 中起主导因素的参数是 𝑇allgather

3.6.3 加速比和效率

加速比：串行并行𝑠(𝑛,𝑝)=𝑇串行(𝑛)𝑇并行(𝑛,𝑝)

效率：串行并行𝐸(𝑛,𝑝)=𝑆(𝑛,𝑝)𝑝=𝑇串行(𝑛)𝑝×𝑇并行(𝑛,𝑝)

p 较小，n 较大的情况下，有近似线性的效率。

3.6.4 可扩展性

• 如果问题的规模能够以一定的速率增加，但效率没有随着进程数的增加而降低，那么该程序就是 可扩展 的。
• 在不增加问题规模的前提下维持恒定效率，此程序称为拥有 强可扩展性。
• 当问题规模增加，通过增大进程数来维持程序效率的，称为 弱扩展性。

3.7 并行排序算法

设总共有 𝑛 个键值，𝑝=comm_sz 个进程，给每个进程分配 𝑛/𝑝 个键值（与前面一样，假定 𝑛 能被 𝑝 整除）。开始时，不对哪些键值分配到哪个进程中加以限制，然而在算法结束时：

• 每个进程上的键值应该以升序的方式存储。
• 若 0≤𝑞<𝑟<𝑝，则分配给进程 𝑞 的每一个键值应该小于等于分配给进程 𝑟 的每一个键值。

所以，如果按照进程编号来进行键值排列（即先是进程 0 的键值，接着进程 1 的，以此类推），所有的键值就可以按升序排列。为了保证表达的清晰性，假设键值都是普通的 int 类型。

3.7.1 简单的串行排序算法

冒泡排序，因为其固有的串行按对排序特性，并行上述算法意义不大。

奇偶排序交换 更适合并行化。关键在于 去耦的比较-交换。算法由一些列阶段组成，这些阶段分 2 种类型：

• 偶数阶段，比较-交换由以下数组对执行： (𝑎[0],𝑎[1]),(𝑎[2],𝑎[3]),(𝑎[4],𝑎[5]),⋯

  而奇数阶段则由以下数组对进行比较-交换： (𝑎[1],𝑎[2]),(𝑎[3],𝑎[4]),(𝑎[5],𝑎[6]),⋯

  这里有个小例子：

  • 开始：5、9、4、3
  • 偶数阶段；比较-交换（5，9）和（4，3），获得序列 5、9、3、4。
  • 奇数阶段；比较-交换（9，3），获得序列 5、3、9、4。
  • 偶数阶段；比较-交换（5，3）和（9，4），获得序列 3、5、4、9。
  • 奇数阶段；比较-交换（5，4），获得序列 3、4、5、9。

  这个例子需要四个阶段来排序四个元素的列表。一般来说，阶段可能会更少些，下面的这个定理保证了我们至多用 𝑛 个阶段排序 𝑛 个元素。

  定理：设 𝐴 是一个拥有 𝑛 个键值的列表，作为奇偶交换排序算法的输入，那么经过 𝑛 个阶段后，𝐴 能够排好序。

void Odd_even_sort(
    int a[]  /* in/out */,
    int n    /* in */
) {
    int phase, i, temp;

    for (phase = 0; phase < n; phase++) {
        if (phase % 2 == 0) { /* Even phase */
            for (i = 1; i < n; i += 2)
                if (a[i - 1] > a[i]) {
                    temp = a[i];
                    a[i] = a[i - 1];
                    a[i - 1] = temp;
                }
        } else { /* Odd phase */
            for (i = 1; i < n - 1; i += 2)
                if (a[i] > a[i + 1]) {
                    temp = a[i];
                    a[i] = a[i + 1];
                    a[i + 1] = temp;
                }
        }
    } /* Odd_even_sort */
}

3.7.2 并行奇偶交换排序

• 任务：在阶段 j 结束时，确定 a[i] 的值
• 通信：确定 a[i] 的值的任务需要与其他确定 a[i+1] 或者 a[i-1] 的任务进行通信，同时，在阶段 j 结束时，a[i] 的值需要用来在阶段 j+1 结束时确定 a[i] 的值。

假设我们有 p = 4 个进程，n = 16 个键值。对分配给每个进程的键值使用快速的串行排序法，如 C 语言库的 qsort 函数对局部链值进行排序。现在，如果每个进程有一个元素，那么进程 0 和 进程 1 交换数据，进程 2 和进程 3 交换数据。我们试着使进程 0 和进程1 交换数据，

Sort local keys;
for(phase = 0; phase < comm_sz; phase++){
    partner = Compute_partner(phase, my_rank);
    if(I'm not idle){
        send my keys to partner;
        receive keys from partner;
        if(my_rank < partner)
            keep smaller keys;
        else
            keep larger keys;
    }
}

Compute_partner：

if (phase % 2 == 0)     /* Even phase */
    if (my_rank % 2 != 0)   /* Odd rank */
        partner = my_rank - 1;
    else                    /* Even rank */
        partner = my_rank + 1;
else                       /* Odd phase */
    if (my_rank % 2 != 0)   /* Odd rank */
        partner = my_rank + 1;
    else                    /* Even rank */
        partner = my_rank - 1;

if (partner == -1 || partner == comm_sz)
    partner = MPI_PROC_NULL;

4 用 Pthreads 进行共享内存编程

从程序员角度，共享内存系统中的任意处理器核都能访问所有的内存区域。因此，协调各个处理器核工作的一个方法，就是把某个内存区域设为 “共享”。

• 对共享内存区域进行更新的代码段，称为 临界区。
• 运行在一个处理器上的一个程序实例称为 线程 (不同于 MPI，在 MPI 中称为进程)。

4.1 进程、线程 和 Pthreads

线程 是轻量级的进程。

• 进程是正在运行（或挂起）的程序的一个实例，除了可执行代码外，它还包括：
  • 栈段。
  • 堆段。
  • 系统为进程分配的资源描述符，如文件描述符等。
  • 安全信息，如进程允许访问的硬件和软件资源。
  • 描述进程状态的信息，如进程是否准备运行或者正在等待某个资源，寄存器中的内容等。
• 线程之间共享一些区域，可以相互访问各自的内存区域。
• 事实上，线程有各自独立的栈和计数器，除此之外，它们基本上共享所有其他区域。

4.2 “Hello, World” 程序

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>

/* 全局变量：所有线程都可获取 */
int thread_count;

void* Hello(void* rank);  /* 线程函数 */

int main(int argc, char* argv[]){
    long       thread;  /* 使用 long 是因为64位系统 */
    pthread_t* thread_handles;
    
    /* 从命令行中获取 线程数量 */
    thread_count = strtol(argv[1], NULL, 10);
    
    thread_handles = malloc(thread_count * sizeof(pthread_t));
    
    for(thread = 0; thread < thread_count; thread++)
        pthread_create(&thread_handles[thread], NULL, Hello, (void*)thread);
    
    
    printf("Hello from the main thread\n");
    
    for(thread = 0; thread < thread_count; thread++)
        pthread_join(thread_handles[thread], NULL);
    
    free(thread_handles);
    return 0;
} /* main */

void* Hell(void* rank){
    long my_rank = (long) rank;
    
    printf("Hello from thread %ld of %d\n", my_rank, thread_count);
    
    return NULL;
} /* Hello */

4.2.1 执行

编译：

gcc -g -Wall -o pth_hello pth_hello.c -lpthread

-lpthread 告诉编译器，我们要连接 Pthreads 线程库。

运行程序：

./pth_hello <number of threads>

./pth_hello 1

4.2.2 准备工作

• 头文件：pthread.h

• strtol 函数功能是将字符串转化为 long int ，它在 stdlib.h 中声明，它的语法形式为：

  C

  折叠 复制

  long strtol(
      const char*  number_p,
      char**       end_p,
      int          base
  );

  • 返回由 number_p 所指向的字符串转换得到的长整数
  • 参数 base 是表达这个整数值所用的基（进制计数制）
  • 如果 end_p 不是 NULL，它就指向 number_p 字符串中第一个无效字符。

4.2.3 启动线程

代码 17 行：

thread_handles = malloc(thread_count * sizeof(pthread_t));

为每个线程的 pthread_t 对象分配内存，pthread_t 数据结构用来存储线程的专有信息。

• pthread_t 对象是一个 不透明 对象，对象中存储的数据都是系统绑定的，用户级代码无法直接访问到里面的数据。

19~21 行代码:

for(thread = 0; thread < thread_count; thread++)
    pthread_create(&thread_handles[thread], NULL, Hello, (void*)thread);

调用 pthread_create 函数来生成线程。它的语法：

int pthread_create(
    pthread_t*             thread_p                 /* out */,
    const pthread_attr_t*  attr_p                   /* in */,
    void*                  (*start_routine)(void*)  /* in */,
    void*                  arg_p                    /* in */
);

• 第一个参数是 pthread_t 对象的指针，必须在调用 pthread_create 函数之前就为 pthread_t 对象分配内存空间。
• 第二个参数不用，所以只是在函数调用时 把 NULL传递给参数。
• 第三个参数表示该线程将要运行的函数。
• 最后一个参数是一个指针，指向传递给 函数 start_routine 的参数。

pthread_create 生成并运行的函数应该是一个类似于下面函数的原型：

void* thread_function(void* args_p);

• 因为类型 void* 可以转换为 C语言任意指针类型，所以 args_p 可以指向一个列表，该列表含有一个或多个 thread_function 函数需要的数值。类似地，thread_function 返回的值也可以是一个包含一个或多个值的列表。

既然线程函数可以接收 void * 类型的参数，我们就可以在 main 函数中为每个线程分配一个 int 类型的整数，并为这些整数赋予不同的数值。当启动线程时，把指向该 int 型参数的指针传递给 pthread_create 函数。

然而，程序员会用类型转换来处理此问题：不是在 main 函数中生成 int 型的进程号，而是把循环变量 thread 转化为 void * 类型，然后在线程函数 hello 中，把这个参数的类型转换为 long 型（第 33 行）。

类型转换的结果是“系统定义”的，但大多数 C 编译器允许这么做。不过，如果指针类型的大小和表示进程编号的整数类型不同，在编译时就会收到警告。在我们使用的机器上，指针类型是 64 位，而 int 型是 32 位，为了避免警告，我们用 long 型替代了 int 型。

4.2.4 运行线程

• 运行 main 函数的线程一般称为主线程。
• 线程的调度是由操作系统来控制的。

4.2.5 停止线程

代码第 25、26行：

for(thread = 0; thread < thread_count; thread++)
    pthread_join(thread_handles[thread], NULL);

为每个线程调用一次 pthread_join 函数，调用一次 pthead_join 将等待 pthread_t 对象所关联的那个线程结束。pthread_join 语法为：