分布式并行计算+MPI

原创已于 2025-05-13 10:30:41 修改 · 1k 阅读

11 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

于 2025-05-10 16:49:21 首次发布

分布式并行计算+MPI

第一章——介绍

一、并行与分布式计算

并行计算：在单一系统中使用两个或多个处理器（计算机）同时工作以解决单个问题
分布式计算：为解决一个计算或者信息处理问题，多台分布的计算机进行协同计算的方式

二、并行与分布式系统的架构

在这里插入图片描述 系统分类：

并行系统（Parallel System）：强调多个处理单元协同完成同一任务。
分布式系统（Distributed System）：强调多个节点通过网络协作。

内存架构：

共享内存（Shared Memory）：所有处理单元可直接访问同一逻辑地址空间，硬件或软件保证一致性。
分布式内存（Distributed Memory）：每个处理节点拥有独立本地内存，必须通过消息传递显式交换数据。

组合	系统分类	内存逻辑模型	编程模型
并行 + 共享	并行系统	共享内存(UMA/NUMA)	OpenMP
并行 + 分布式	并行系统	分布式内存	MPI
分布式 + 分布式	分布式系统	分布式内存	消息队列
分布式 + (虚拟)共享	分布式系统	虚拟共享内存	DSM 库

三、并行计算的必要性

技术驱动的必然

摩尔定律转向：单核频率提升受限，现在更多依靠增加核心数量来提升性能。
功耗与能效约束：提高频率会大幅增加能耗，而多核可以在低频下同时工作，实现更高能效。
指令级并行（ILP）极限：传统靠硬件自动提高程序运行效率的方法已接近上限。

应用需求的驱动

科学仿真
人工智能

四、并行编程

并行编程是指通过编程语言显式地指示程序中不同部分在多个处理器上并发执行。请添加图片描述

常见的并行化方法

重写并行程序：直接编写能够在多核系统上运行的代码。
编译指示技术：如 OpenMP 提供的 #pragma，指导编译器生成并行代码。
自动并行化：依靠编译器分析串行代码，自动转换为并行形式。

并行模式

数据并行：每个核心使用相同的操作处理不同的数据子集
任务并行：不同核心执行不同的任务

并行编程的挑战

同步：确保多个核心协调一致工作，不出错
通信：核心间需要传递中间结果
负载均衡：确保每个核心处理的任务量大致相同，避免资源浪费

第二章 ——并行体系结构

一、弗林分类法

类别	含义	示例
SISD	单指令流单数据流（传统顺序处理器）	单核 CPU
SIMD	单指令流多数据流	GPU、向量处理器
MISD	多指令流单数据流（理论性强）	极少见，用于容错处理
MIMD	多指令流多数据流	多核处理器、分布式系统、大多数并行机

弗林分类法关注处理器执行模型，本质是从处理器执行逻辑层面出发，分析系统在时间某一时刻能同时处理多少条指令/数据。

编程模型	弗林分类	系统结构层级说明
单核 CPU	SISD	单处理器 + 共享内存
CUDA	SIMD	并行系统 + 共享内存（NUMA）
OpenMP	MIMD	并行系统 + 共享内存（UMA/NUMA）
MPI	MIMD	并行系统 + 分布式内存

SIMD(单指令流、多数据流)

在这里插入图片描述

有一个控制单元发出统一指令，有多个处理单元分别执行该指令但作用于不同数据，常见架构有向量处理器、GPU等。

优点：

在数据并行场景下高效
适合结构化数据处理

局限性：

不适合任务并行场景
不支持复杂控制逻辑和数据依赖
“分支发散”问题

MIMD(多指令流、多数据流)

在这里插入图片描述

每个核心都有独立的控制单元和数据路径支持异步执行，可使用共享内存或分布式内存。

优点：

适应任务并行、数据并行等多种场景
可扩展性强
支持复杂控制逻辑和数据依赖

局限性：

编程复杂
通信开销大
对缓存一致性和负载均衡要求高

二、共享/分布式内存并行系统

共享内存系统

在这里插入图片描述

多个处理器通过互连网络连接到统一的存储系统；所有处理器**可以直接访问所有内存数据。

UMA

在这里插入图片描述

NUMA

在这里插入图片描述

类型	UMA（均匀存储访问）	NUMA（非均匀存储访问）
存储结构	所有处理器通过互联网络访问统一主内存	每个处理器连接一个本地内存块，形成分布式共享内存
访问延迟	均匀：所有处理器访问任何地址延迟相同	非均匀：访问本地内存快，访问远程内存慢
缓存系统	每个处理器可有私有缓存	每个处理器同样带私有缓存

分布式内存系统

在这里插入图片描述

每个处理器具有独立的私有内存，没有全局统一地址空间；处理器之间不能直接访问对方的内存；通信必须通过显式的消息传递机制。

区别与讨论

维度	共享内存系统	分布式内存系统
地址空间	统一地址空间	分离地址空间
通信方式	隐式（通过共享变量）	显式（通过消息传递）
可扩展性	较差（通信总线成为瓶颈）	较强（节点之间独立，易于扩展）

共享内存系统与分布式内存系统的本质区别是是否为统一地址空间，因此虽然我们看NUMA与分布式内存系统很像，是因为它们在一样在物理上内存是分开的，但是NUMA将分布的内存看作一整个内存从而统一地址空间。因为共享内存系统每个CPU都能看到一整个内存，加上它们可以有自己的私有缓存，因此就有了缓存一致性的问题。分布式内存系统，本节点的内存不会被其它节点直接修改，需要消息传递通知该节点修改，同时该节点修改也会消息传递修改，因此缓存一致。

三、互联网络——Interconnect

互联网络是并行系统中各处理器与内存模块等之间数据传输机制。其实考虑到并行系统存在共享内存系统和分布式内存系统，这里互联网络不止包括计算机内部的互联还包括计算机之间的网络。

集中式互联（共享内存互联）

类型	描述	特点	典型应用
总线	所有处理器通过一条共享总线访问内存	成本低、冲突高、扩展性差	小型SMP、多核芯片
交叉开关	每个处理器与内存模块之间有专用通路	并发高、成本高、连接复杂	高性能SMP、多插槽服务器

分布式互联（分布式内存互联）

类型	描述	特点
静态/直接互联	拓扑结构固定，通信路径不随负载变化	路由简单、易建模，适合规则型并行结构
动态/间接互联	利用交换器动态建立通信路径，适应性更强	适合大规模异构通信、灵活、调度复杂

值得注意的是上面的分布式互联也用于NUMA共享内存系统。

静态/直接互联

指标

节点度：每个节点直接连接的邻居节点数
网络直径：任意两节点之间的最长最短路径长度
对分宽度：将网络划分为两半，最少需切断的边数
对分带宽：对剖宽度 × 单条链路带宽（即实际吞吐能力）
是否恒定边长：网络中任意两相连节点之间的物理链路距离是否一致

全连接网络

每个节点都直接连接到其他每个节点。全连接成本太高，实际实现不可能，但可以作为理论边界。

节点度： $\in O(N)$
网络直径： $\in O(1)$
对分宽度： $N24∈O(N2)\frac{N^{2}}{4} \in O(N^2)$
是否恒定边长：否，物理上，为了连接每对节点，就必须铺设一条物理链路，如果将所有节点排列在二维平面或三维空间中，一些节点之间的物理距离会远远大于另一些。

线状网络

交换机排列成一个一维网状。

节点度： $\in O(1)$
网络直径： $\in O(N)$
对分宽度： $\in O(1)$
是否恒定边长：是，物理上可以保证两个节点之间链路一样长。

环状网络

线状网络变体，交换机连接成环状。

节点度： $\in O(1)$
网络直径： $⌊N2⌋∈O(N)\lfloor\frac{N}{2}\rfloor \in O(N)$
对分宽度： $\in O(1)$
是否恒定边长：是，物理上可以将节点均匀排列在圆周上，每条边等长。

2D Mesh 网格

在这里插入图片描述

节点排列成一个二维的格子或网格，只允许在相邻的交换机之间进行通信。

节点度： $\in O(1)$
网络直径： $2(N−1)∈O(N12)2(\sqrt{N} - 1) \in O(N^\frac{1}{2})$
对分宽度： $N∈O(N12)\sqrt{N} \in O(N^\frac{1}{2})$
是否恒定边长：是，因为每个节点只连接到邻接节点，这些链路在物理布局中长度一致，易于实现布线。

2D Torus 环面

在这里插入图片描述

2D Mesh 网格变体，包含了在网格边缘之间的环绕连接。

节点度： $\in O(1)$
网络直径： $2⌊N2⌋∈O(N12)2\lfloor\frac{\sqrt{N}}{2}\rfloor \in O(N^\frac{1}{2})$
对分宽度： $2N∈O(N12)2\sqrt{N} \in O(N^\frac{1}{2})$
是否恒定边长：否，可以若是一维环形尚且可以做成圆柱实现物理链路相等，二维环面难以实现链路相等。

超立方体网络

在这里插入图片描述

超立方体节点编号为 $d$ 位二进制数（共有 $2^d$ 个节点），每个节点与仅有一位不同的其他 $d$ 个节点相连，本质上是 $d$ 个一位不同的邻居，每一维提供一条边。

节点度： $\in O(logN)$
网络直径： $\in O(logN)$
对分宽度： $\in O(N)$
是否恒定边长：否，因为高维嵌入时，连接的节点在物理距离上可能差异很大；例如 8 维超立方体（256 个节点）中，一维连接可能对应物理上很远的两个节点。

树状网络

几种特殊的二叉树 - 满二叉树、完全二叉树、二叉排序树、平衡二叉树 - csy123 - 博客园

树形网络通常为完全二叉树或 $k$ -叉树， $N=∑i=1hki−1=kh−1k−1N=\sum^{h}_{i=1}k^{i-1}=\frac{k^{h}-1}{k-1}$

节点度： $\in O(logN)$
网络直径： $2log_{k}N \in O(logN)$
对分宽度： $\in O(1)$
是否恒定边长：是，可以看作有分支的线状网络，因此每个节点之间的物理链路可以做到相等。

动态/间接互联

在这里插入图片描述

消息传输时间

在这里插入图片描述

时延：从源头开始传输数据到目的地开始接收第一个字节之间所经过的时间
带宽：目的地在开始接收第一个字节后接收数据的速率

蝶形网络

在这里插入图片描述

蝶形网络就像是超立方体的展开版本，将超立方体的按位变化路径具体展开。若采用d位编码，则存在 $2^{d}$ 个处理器和内存节点， $d*2^{d}$ 个交换节点。

节点 $(i, x)$ 向下一层 $(i + 1, x)$ 和 $(i + 1, x^{'})$ 连接， $x^{'}$ 是 $x$ 的第 $i$ 位（从左到右）翻转后的结果，即每一层决定翻转某一位。

树状网络

在动态树状网络中，中间层的节点是交换节点，叶子节点是处理器和内存。

Fat Tree Network

在这里插入图片描述

四、缓存一致性

概念

缓存一致性确保多个处理器各自缓存中对同一内存地址的副本始终保持一致性。否则，程序的并行执行将出现不符合预期的行为。

虚假共享(false sharing)，当两个不同处理器访问的是不同变量，但这些变量恰好共享同一个缓存行，只要一个处理器写入，另一个就会因一致性协议失效其缓存行。虚假共享并不会导致不正确的结果，但可能会破坏一个程序的性能，因为它导致对内存的访问比预期的多得多。

缓存写回策略	说明
写回（Write-back）	修改只写入缓存，等缓存替换时或被失效时才写入主内存
写直达（Write-through）	每次写操作立即同时写入缓存和主内存

缓存写回策略描述缓存与主内存之间的数据同步时机，目标是保证主内存与缓存之间的数据一致性。

缓存一致性策略	说明
写失效(Write Invalidate)	写操作先让其他核的对应缓存行失效，自己变为唯一写者
写更新(Write Update)	写操作时直接将新值广播给其他缓存，让它们更新副本
目录机制(Directory-based)	每块数据都有一个目录记录拥有副本的处理器

缓存一致性策略解决多核缓存副本之间的一致性，目标是确保其他处理器不使用旧值。

现代系统更常用Write Invalidate，因为写操作通常是少数，同时保持写者简单，不需广播内容。

常见一致性协议

MSI 协议

三个状态：Modified（已修改）/ Shared（共享）/ Invalid（无效）；
是最基础的写无效协议，读写引发状态切换。
更适合UMA；
有进一步发展的MESI(增加Exclusive状态，实际中更常用)、MOESI(增加Owned状态)

Directory-based 协议

中心目录跟踪每个块的状态和所有缓存副本；
更适合NUMA ；
避免广播，减少总线开销。

基于监听的缓存一致性

基本原理：所有缓存通过共享总线监听（snoop）通信，任何一个处理器发出对内存的读写请求时，其他缓存都能感知该操作并响应（失效/更新等）。

基于监听的缓存一致性使用于共享总线结构（如 UMA），一般采用写失效协议（如 MSI、MESI、MOESI）。

工作流程（写回 + 写失效）

核 A 对地址 X 写命中，并广播；
其他缓存监听总线，若持有 X 标为Invalid；
A 缓存将 X 标为 Modified；
写操作只修改缓存，主存暂不更新；
当缓存行被替换或失效时，才写回主存。

基于目录的缓存一致性

状态	含义
U（Uncached）	该块未缓存在任何处理器中。目录中没有标记此块对应的处理器位。
S（Shared）	该块被一个或多个处理器以只读方式缓存。目录中记录所有拥有副本的处理器编号（位图）。
E（Exclusive）	该块被一个处理器以独占方式缓存并可能已修改。

触发事件	当前状态	动作	新状态
读请求	`U`	内存提供数据；记录请求者	`S`
读请求	`S`	内存提供数据；更新 Presence Vector	`S`
读请求	`E`	请求独占者写回并降级为共享；发副本	`S`
写请求	`S`	向所有共享副本发送 Invalidate；设新拥有者	`E`
写请求	`E`	如果拥有者≠请求者，则转移独占权限；或直接写	`E`

第三章——并行程序设计

一、Foster 设计方法论

阶段	关注问题	是否与机器相关	目标
划分	任务/数据独立性	否	最大化并发性
通信	数据交换依赖	否	最小化通信需求
聚合	并行单元组合	是	提高任务粒度，降低通信频率
映射	任务到处理器绑定	是	优化实际执行性能

前两步是在“纸面上设计”，决定程序逻辑层面的并行结构。后两步是“具体部署”，决定如何在硬件平台上高效执行程序。

划分

划分将问题分解成大量的基本任务（子问题），是并行设计的核心第一步，目标是识别出可以并行执行的最小计算单元。策略有域分解（数据划分）和功能分解（任务划分），这也对应着第一章中提到的数据并行和任务并行，这两者不是互斥的而是互补的。

域分解，举个例子，假设有两个N阶矩阵相乘，我们可以将第一个矩阵换分为行的数组，将第二个矩阵划分为列的数组，这样每个节点只要计算一行一列的向量乘法即可。功能分解，举个例子，假设有一个公式要求先将数进行平方，再将这个数哈希，我们可以有两个节点，第一个节点负责平方，第二个节点负责哈希，这样当第一波数据进行哈希时第二波数据可以进行平方，通过流水线进行并行处理。

通信

通信是明确划分后任务之间需要共享或交换的数据，并设计通信方式，目标是确定任务之间的数据依赖关系。类型有本地通信和全局通信。

为什么对于功能分解来说，明确通信相对容易，但对于域分解来说却很困难？
功能分解后功能块清晰，数据流显式，谁产生谁消费很明确；域分解后数据分布在多个任务中，依赖隐含在空间结构中（如边界数据共享）。

本地通信数据交换只发生在邻近任务之间，通常是结构化的、规律性的通信，例如2D 网格中每个任务与上下左右邻居交换边界数据。
全局通信涉及任务集合中的大范围数据交互，通常不在早期设计时显式建通道，缺点是通信量大、易造成同步瓶颈。

全局通信策略示例（优化方式）

集中式求和每个任务把一个值发给根任务（S），S 汇总，简单但不可并行，且 通信瓶颈严重。
流水线式前缀和将 N 个任务线性排列为数组，每个任务 i 接收 i-1 的部分和，加上自己的值，多次执行可流水线并行。
树形结构的分治求和问题被递归划分为子问题，结果通过树形结构归并，具有规则性、可并行性强，是结构化的全局通信策略。

为什么流水线和分治法不是“本地通信”？
本地通信要求数据只在空间局部范围流动；全局通信要求涉及所有（或大多数）任务的结果聚合或广播。流水线式和分治法式结构化地汇聚所有结果，属于全局通信。

聚合

聚合将划分和通信中得到的基本任务和通信组合成更大、更高效的复合任务。在并行系统中任务过小会导致调度/通信成本变高因此性能受限，聚合目的是控制并行粒度，降低调度和通信成本。

聚合的方式

在这里插入图片描述

合并通信双方任务（如图 a），如果任务 A 输出数据给 B，且 A 和 B 执行时间顺序固定，可将 A 与 B 合并为一个复合任务使通信内部化，通信代价消失。
合并多个通信方向的任务为组（如图 b），将多个任务（发送者/接收者）合并成“组”，共享通信结构减少通信数量和调度次数。

矩阵案例

操作一个大小为 8×128×256 的三维矩阵在不同规模的 集中式多处理器（SMP） 上运行

Q1：将第2和第3维聚合（即任务变成 8 个 128×256 区块） ➜ 能否在4核机器上运行？
生成 8 个任务，可映射到 4 核机器（每核运行 2 个任务），聚合减少了任务数量，通信更局部

Q2：如果目标机器有 8 个 CPU？
8 个任务刚好映射，每个任务独占一个 CPU，最优分配

Q3：如果目标机器有超过 8 个 CPU？
任务数固定为 8，额外CPU闲置，需要更细粒度的划分（取消聚合或更细拆）

映射

映射是将已经设计好的并行任务，分配到实际硬件资源（处理器/核心）上的过程。

在这里插入图片描述

目标	描述
最大化处理器利用率	保证每个处理器的工作负载均衡，避免闲置（= 负载均衡 / 计算平衡）
最小化通信代价	把通信频繁的任务映射到物理上靠近的处理器（或同一处理器），减少通信成本

两者存在矛盾，完全平衡负载，可能导致跨节点频繁通信，强调通信局部性，可能让某些核任务过重。因此映射是一个NP问题，必须依赖于启发式算法找到相对最优的方案。

映射的基本策略有静态映射和动态映射。静态映射编译时就分配好任务到处理器，常见于结构规则、负载已知的问题；优点是效率高，开销低，缺点是无法适应负载变化。动态映射运行时分配任务，适合负载不均或不确定的问题；优点是可适应性好，缺点是调度开销更大，通信更复杂。

讨论

前面说，前两步是在“纸面上设计”，决定程序逻辑层面的并行结构，后两步是“具体部署”，决定如何在硬件平台上高效执行程序。实际上第三步聚合也是逻辑层面的问题，但是它“间接考虑”物理硬件（用作粒度设计的依据），所以也可以认为它是硬件平台的问题。可以这样说，划分和通信是找到理论中基本并行任务与对应的通信关系，聚合是将理论结合到实际，调整任务粒度，最后映射将其具体实施到硬件上。

阶段	目标	输出
划分	把问题拆成子任务	基本任务
通信	明确任务间数据依赖	通信关系
聚合	合并任务减少通信，调整任务粒度	一组复合任务和更简洁的通信结构
映射	把复合任务分配给处理器	最终的“任务–资源”调度计划

二、案例研究

1.一维热传导问题（边值问题或扩散问题）

在这里插入图片描述

上面两幅图是这个问题的直观了解，一根金属杆长一个单位，它的温度只受两端的冰水的影响，初始状态是两端为0摄氏度，中间为100摄氏度的抛物线。

在这里插入图片描述

并行程序设计

划分：划分是把问题拆成子任务，一整根金属杆可以看作N个金属段，温度的连续变化可以看作dt时间均匀变化的累计结果，因此问题被划分为N个金属段的dt时间的变化，如上图（a）所示。
通信：一个金属段的温度变化只受两端金属段和自身温度的影响，因此每dt时间迭代需要来自两端的温度和自身温度的信息，如上图（b）所示每个金属段有两个接收通信和两个发送通信（自身不需要通信），这是典型的本地通信的例子。
聚合：如上图（c）所示，将10个金属段节点聚合为4：3：3的3个复合任务（这里考虑机器只有3个处理器），因此将每次迭代需要的18次通信缩减为4次通信，并且恰好可以每个处理器只处理一个复合任务，避免频繁调度。
映射：这里将3个复合任务分配到对应的3个处理器上即可。

时间复杂度分析

假设单位金属杆分为N段，时间间隔dt，迭代K次

迭代公式 $uit+1=uit+αui−1t−2uit+ui+1t(Δx)2Δtu^{t+1}_{i}=u^{t}_{i}+\alpha\frac{u^{t}_{i-1}-2u^{t}_{i}+u^{t}_{i+1}}{(\Delta x)^{2}}\Delta t$
显然每次迭代可以常数步计算出来，因此迭代部分代码时间复杂度 $O (1)$

串行算法每次迭代需要更新N-2个值(假设端点值恒为0)，因此迭代一次时间复杂度为 $O (N)$ ，那迭代K次时间复杂度为 $O (N K)$ 。

并行算法每次迭代每个节点需要更新约 $NP\frac{N}{P}$ 个值，因此迭代一次时间复杂度 $O(NP)O(\frac{N}{P})$ ，假设通信时间复杂度为 $O(N^{'})$ ，那么K次迭代时间复杂度为 $O(K(NP+N′))=O(NKP+N′K)O(K(\frac{N}{P}+N^{'}))=O(\frac{NK}{P}+N^{'}K)$ 。

如果要计算执行时间，假设 $μ\mu$ 为更新值时间， $λ\lambda$ 为通信时间

$Tserial=K(N−2)μT_{serial}=K(N-2)\mu$ ； $Tparallel=K(NPμ+2λ)T_{parallel}=K(\frac{N}{P}\mu+2λ)$

2.规约问题

规约（Reduction）问题是并行计算中一种常见的操作，主要指将一组数据元素使用某个结合操作符 $⊕\oplus$ 聚合成一个单一值的过程，只要操作符 $⊕\oplus$ 满足结合律 $\oplus b) \oplus c = a \oplus (b \oplus c)$ ，就可以用并行方式高效计算。

操作符（⊕）	含义	示例
+	求和	a₀ + a₁ + … + aₙ₋₁
×	求积	a₀ × a₁ × … × aₙ₋₁
∧	位与	a₀ ∧ a₁ ∧ … ∧ aₙ₋₁
∨	位或	a₀ ∨ a₁ ∨ … ∨ aₙ₋₁
max	求最大值	max(a₀, a₁, …, aₙ₋₁)
min	求最小值	min(a₀, a₁, …, aₙ₋₁)

并行程序设计

在这里插入图片描述

假设我们有 $N=2^{d}$ 个数据进行结合计算。

划分：将N个数据的结合计算根据结合律不断加括号，最后每个括号里面只有两个值进行结合计算，如 $(x0⊕x1)⊕(x2⊕x3)(x_0 \oplus x_1) \oplus (x_2 \oplus x_3)$ ，形成对数深度任务树，这样共有 $2^{d-1}$ 个基本任务，其实也可以是 $2^{d}$ 个基本任务，但这样第一轮没有直接通信，为了便于理解增加通信无意义，因此采用 $2^{d-1}$ 。
通信：为 $2^{d-1}$ 个基本任务使用 $d - 1$ 位编码，这样可以使用超立方体网络进行通信，进行 $d - 1$ 次迭代，第 $i$ 次可以使第 $i$ 位为0的对为1的通信计算结果，为1的收到通信最后将通信数据再与自身计算后的结果结合计算，经过所有迭代后，在编码全为1的节点就可以得到最终计算结果，属于典型的全局通信。
聚合：假设d为4，处理器数量p为4，如上图所示，可以将8个基本任务聚合为4个复合任务，复合任务仍可以采用超立方体网络进行通信。
映射：将4个复合任务分配到4个处理器上。

时间复杂度分析

$Tserial=(N−1)μ=(2d−1)μ∈O(d)T_{serial}=(N-1)\mu=(2^{d}-1)\mu\in O(d)$

$Tparallel=(NP−1)μ+(logP−1)(μ+2λ)∈O(NP+logP)T_{parallel}=(\frac{N}{P}-1)\mu+(logP-1)(\mu+2\lambda) \in O(\frac{N}{P}+logP)$

基本算法

基本算法的两类抽象：Scatter 与 Gather

类别	方向	功能	通信结构
Gather	多 → 一 / 多	聚合数据、收集结果	树 / 超立方体
Scatter	一 → 多	分发数据、广播初始化	树 / 超立方体

Scatter与Gather两类操作在拓扑结构中对称，因此常常有相同的优化通信方法。

Gather 的结构化形式：

Gather：树状结构从底向上收集；
All-Gather：每轮互发互收，常在环/超立方体中实现。
Reduce：Gather + 结合计算（例如 $⊕\oplus$ ），如规约；
All-Reduce：结构化 All-Gather + Reduce，如规约双向通信规约。

Scatter 的结构化形式

Scatter：树状结构从顶向下分发；
Broadcast：一种特殊的 Scatter（所有进程都接收同一份数据）；
结构化 Scatter：每个节点仅传自己负责的那一部分数据给下层；

完全图（全连接网络）vs 超立方体网路

假设有数据总量 $n$ ，节点 $p$ ，通信时间 $λ\lambda$ ，带宽 $β\beta$ ，每个节点保存 $np\frac{n}{p}$ 的数据，要求是的一个节点获得全量数据（Gather）。

$Tfull=(p−1)(λ+npβ)T_{full}=(p-1)(\lambda+\frac{n}{p\beta})$

$Thypercube=∑i=1logp(λ+2i−1npβ)T_{hypercube}=\sum^{logp}_{i=1}(\lambda+\frac{2^{i-1}n}{p\beta})$

为什么超立方更优?

比较维度	完全图	超立方体	哪个更优
通信轮数	$p - 1$ （线性）	$log⁡p\log p$ （对数）	超立方
通信端口数	每节点连接 $p - 1$ 节点	每节点连接 $log⁡p\log p$ 个	超立方
扩展性	不易扩展（连接数爆炸）	易于扩展	超立方
总通信量	固定，线性增长	分层增长，整体低一些	超立方（在大p时）

第四章——性能评估

一、加速比和效率

加速比及理论上界

$\frac{Ts}{Tp}$

这个公式加速比计算公式，其中 $T s$ 为串行算法运行时间， $T p$ 为使用 $p$ 个处理器并行算法运行时间。

在这里插入图片描述

类型	数学关系
线性加速比	$S_p = p$
超线性加速比	$S_p > p$
亚线性加速比	$S_p < p$

超线性加速比不是错误，而是说明串行程序未充分利用硬件（如缓存）。但从理论上，线性加速比已是最理想状态，超线性通常表示性能上的“意外红利”。

$ψ(n,p)≤σ(n)+ϕ(n)σ(n)+ϕ(n)/p+κ(n,p)\psi(n,p)\leq \frac{\sigma(n)+\phi(n)}{\sigma(n)+\phi(n)/p+\kappa(n,p)}$

符号	含义
$ψ(n,p)\psi(n, p)$	加速比的实际上限，问题规模为 $n$ ，处理器数为 $p$
$σ(n)\sigma(n)$	串行计算时间（无法并行的部分）
$ϕ(n)\phi(n)$	完全可并行的计算部分串行计算时间
$κ(n,p)\kappa(n, p)$	并行开销

这个公式是加速比理论上界表达式，它比经典的 $Sp=TsTpS_p = \frac{T_s}{T_p}$ 更精细地建模了实际并行程序的开销，适用于考虑实际通信、同步和串行负载等开销的复杂模型。

当处理器数量 $p$ 增加时，加速比通常会经历先快速上升、后趋于平缓甚至下降的阶段。这通常会出现一个拐点。在拐点增加处理器带来的性能提升开始无法弥补额外的开销，是设计系统时计算资源与性能收益的平衡点。

效率及理论上界

$Ep=SppEp=\frac{Sp}{p}$

效率表示并行程序实际利用了多少处理器能力，也就是每个处理器贡献的有效计算比例。

加速比类型	效率	说明
线性加速比	$E_p = 1$	100% 效率，理想情况
亚线性加速比	$E_p < 1$	常见，存在开销或负载不均
超线性加速比	$E_p > 1$	稀有，表示额外性能收益

$\epsilon(n,p) \leq \frac{\sigma(n)+\phi(n)}{p\sigma(n)+\phi(n)+p\kappa(n,p)} $ $\leq \epsilon(n,p) \leq 1$

最小值0为极端低效，所有时间花在开销上，或任务过小；最大值1为理想情况，无串行部分，也无开销，即线性加速。

二、指标

Amdahl定律

$ψ(n,p)≤σ(n)+ϕ(n)σ(n)+ϕ(n)/p+κ(n,p)≤σ(n)+ϕ(n)σ(n)+ϕ(n)/p\psi(n,p)\leq \frac{\sigma(n)+\phi(n)}{\sigma(n)+\phi(n)/p+\kappa(n,p)} \leq \frac{\sigma(n)+\phi(n)}{\sigma(n)+\phi(n)/p}$

Amdahl定律从加速比的理论上界出发，忽略通信、调度等并行开销，假设问题规模固定，比如处理100张图片，只通过增加处理器数量来提高速度，用于估算最大可能加速比，即使用更多处理器能快多少倍完成任务。

令 $\frac{\sigma(n)}{\sigma(n)+\phi(n)}$ ，那么 $ψ(n,p)≤1f+(1−f)/p\psi(n,p)\leq \frac{1}{f+(1-f)/p}$ ，这里 $f$ 是程序中无法并行的串行比例 $\leq f \leq 1$ 。

可以观察到的是串行部分 $f$ 限制了加速比的上限，即使你有再多处理器，只要有一小部分是串行的，就不可能无限加速。最大可能加速比为 $lim⁡p→∞Sp=1f\lim_{p \to \infty} Sp = \frac{1}{f}$ 。

例题

程序95%的执行时间发生在一个可以并行执行的循环内。如果使用8个CPU执行程序的并行版本，期望的最大加速比为多少？
$Spmax=10.05+0.95/8≈5.9Sp_{max}=\frac{1}{0.05+0.95/8}\approx5.9$

Gustafson-Barsi定律

$ψ(n,p)≤σ(n)+ϕ(n)σ(n)+ϕ(n)/p+κ(n,p)≤σ(n)+ϕ(n)σ(n)+ϕ(n)/p\psi(n,p)\leq \frac{\sigma(n)+\phi(n)}{\sigma(n)+\phi(n)/p+\kappa(n,p)} \leq \frac{\sigma(n)+\phi(n)}{\sigma(n)+\phi(n)/p}$

Gustafson-Barsi定律同样从加速比的理论上界出发，忽略通信、调度等并行开销，假设运行时间固定，比如现在有1个小时可以处理图片，通过增加处理器数量来完成更多工作，用于估算最大可能比例加速度，即一个小时里面使用更多处理器能处理多少倍的图片。为什么从加速比的公式出发呢，将 $σ(n)\sigma(n)$ 看作无法并行的串行部分运行时间，那么剩余时间p处理器可以并行处理的问题规模在单处理器上运行时间为 $ϕ(n)\phi(n)$ ，而单处理器串行处理的问题规模运行时间 $ϕ(n)/p\phi(n)/p$ ，这里分子分母使用在单处理器上串行运行的时间量化问题规模，因此可以使用加速比的公式来定义比例加速度。

令 $s=σ(n)σ(n)+ϕ(n)/ps=\frac{\sigma(n)}{\sigma(n)+\phi(n)/p}$ ，那么 $ψ(n,p)≤p+(1−p)s\psi(n,p)\leq p+(1-p)s$ ，这里 $s$ 表示在固定时间内，并行程序处理串行部分运行时间占总时间的比例，因为假设运行时间固定，因此并行程序和串行程序总时间为 $σ(n)+ϕ(n)/p\sigma(n)+\phi(n)/p$ ，因此并行程序处理串行部分运行时间占总时间的比例为 $s$ 。

例题1

一个运行在10个处理器上的应用程序在串行代码中花费了3%的时间。该应用程序的比例加速度是多少？
$ψ=10+(1−10)∗0.03=10−0.27=9.73\psi=10+(1-10)*0.03=10-0.27=9.73$

例题2

如果一个程序要在8个处理器上达到7的比例加速度，那么并行程序在串行代码中可以花费的最大比例是多少？
$7 = 8 + (1 - 8) s$ $s=17≈0.14s=\frac{1}{7}\approx0.14$

Karp-Flatt指标

Amdahl 定律和 Gustafson-Barsis 定律都忽略了通信、调度等并行开销，也就是 $κ(n,p)\kappa(n, p)$ 这一项的影响。在实际系统中，这些开销会导致加速比（Amdahl）或比例加速度（Gustafson）被高估，不能准确反映真实的并行性能。

为了更准确评估程序中不可并行的瓶颈部分，Karp 和 Flatt 提出了另一个度量指标：Karp–Flatt 指标 $e (p)$ ，它通过实际测得的加速比反推程序的“串行性”，能更好地考虑总系统开销。

$e=σ(n)+κ(n,p)σ(n)+ϕ(n)=1ψ−1p1−1pe=\frac{\sigma(n)+\kappa(n,p)}{\sigma(n)+\phi(n)}=\frac{\frac{1}{\psi}-\frac{1}{p}}{1-\frac{1}{p}}$

这个指标使用了加速比的上界公式推导，这里要打个补丁，其实这个上界公式理论上就是n、p下的加速比，因此可以代入指标推导，不过与实际有些误差，但是仍可以通过这个指标反应程序的并行性。

例题1

$p$	2	3	4	5	6	7	8
$ψ\psi$	1.8	2.5	3.1	3.6	4.0	4.4	4.7

在8个CPU上加速比仅为4.7的主要原因是什么？

$e$	$19\frac{1}{9}$	$110\frac{1}{10}$	$331\frac{3}{31}$	$772\frac{7}{72}$	$110\frac{1}{10}$	$13132\frac{13}{132}$	$33329\frac{33}{329}$

$e≈0.1e\approx0.1$ ，由于常数e是固定的，因此较大的串行部分是导致加速比较低的主要原因。

例题2

$p$	2	3	4	5	6	7	8
$ψ\psi$	1.9	2.6	3.2	3.7	4.1	4.5	4.7

在8个CPU上加速比仅为4.7的主要原因是什么？

$e$	$119\frac{1}{19}$	$113\frac{1}{13}$	$112\frac{1}{12}$	$13148\frac{13}{148}$	$19205\frac{19}{205}$	$554\frac{5}{54}$	$33329\frac{33}{329}$
$≈\approx$	0.052	0.076	0.083	0.878	0.927	0.926	0.100

由于常数e不断增加，因此并行开销是导致加速比较低的主要原因。

一些结论

当 $e (p)$ 下降：程序具备良好并行性，扩展性强。
当 $e (p)$ 接近常数：程序的串行部分是主要瓶颈。
当 $e (p)$ 增加：并行开销显著，系统瓶颈来自通信或调度等因素。

例题3

$p$	4	8	12
$ψ\psi$	3.9	6.5	？

这个程序在 12 个处理器上是否有可能达到 10 倍加速比？

$e$	$1117\frac{1}{117}$	$391\frac{3}{91}$

$e (p)$ 显著上升，说明并行开销随着处理器数量增加变得更为显著，继续增加处理器会进一步增加通信与同步开销，导致加速比增长趋缓，我们可以估算 12 个处理器下的理论加速比上限：
$\psi(12) = \frac{1}{e(8)(1 - \frac{1}{12}) + \frac{1}{12}} = \frac{1}{0.033 \cdot \frac{11}{12} + \frac{1}{12}} \approx 7.5$
因此，该程序在 12 个处理器上不可能达到 10 倍加速比

可扩展性

并行系统的可扩展性是衡量其随着处理器数量增加而提高性能的能力。可扩展的系统可以保持效率，即在增加处理器时仍然能够保持高效率。

等效率性是衡量可扩展性的一种方式。

推导等效关系

并行开销公式： $To(n,p)=(p−1)σ(n)+pκ(n,p)T_o(n,p)=(p−1)σ(n)+p\kappa(n,p)$ ，指的是所有进程花在原有串行算法以外操作的全部时间。

将并行开销代入加速比公式： $ψ(n,p)≤σ(n)+ϕ(n)σ(n)+ϕ(n)/p+κ(n,p)≤p(σ(n)+ϕ(n))σ(n)+ϕ(n)+To(n,p)\psi(n,p)\leq \frac{\sigma(n)+\phi(n)}{\sigma(n)+\phi(n)/p+\kappa(n,p)} \leq \frac{p(\sigma(n)+\phi(n))}{\sigma(n)+\phi(n)+T_o(n,p)}$

得到效率公式： $ϵ(n,p)=ψp≤σ(n)+ϕ(n)σ(n)+ϕ(n)+To(n,p)\epsilon(n,p)=\frac{\psi}{p}\leq\frac{\sigma(n)+\phi(n)}{\sigma(n)+\phi(n)+T_o(n,p)}$

由于 $T(n,1)=σ(n)+ϕ(n)T(n,1)=σ(n)+\phi(n)$ ，于是 $ϵ(n,p)≤11+To(n,p)T(n,1)\epsilon(n,p)\leq\frac{1}{1+\frac{T_o(n,p)}{T(n,1)}}$

当效率不变时， $E≤11+To(n,p)T(n,1)E\leq\frac{1}{1+\frac{T_o(n,p)}{T(n,1)}}$

因此 $T(n,1)≥E1−ETo(n,p)T(n,1)\geq\frac{E}{1-E}T_o(n,p)$ ，即 $T(n,1)≥CTo(n,p)T(n,1)\geq CT_o(n,p)$ ，这就是等效关系

当问题规模n增大时，必须保证串行运行时间 $T (n, 1)$ 的增长速度不低于并行开销 $T_o(n,p)$ ，否则并行效率会下降。

可扩展性函数

上面得出等效关系 $T(n,1)≥CTo(n,p)T(n,1)\geq CT_o(n,p)$ ，由于n 增大时，必须保证串行运行时间 $T (n, 1)$ 的增长速度不低于并行开销 $T_o(n,p)$ ，不妨将其化简为 $\geq f(p)$

令 $M (n)$ 问题规模为n时所需的内存， $M (f (p)) p$ 表明每个处理器的内存用量必须增加才能保持相同的效率，称 $M (f (p)) p$ 为可扩展性函数。

为了在增加p时保持效率不变，我们必须增加n，最大的问题规模受限于可用内存，它与p是线性关系，可扩展性函数显示了每个处理器的内存使用量必须如何增长以保持效率不变，可扩展性函数为常数意味着并行系统是完美可扩展的。

在这里插入图片描述

例题1——规约

串行算法复杂度为 $\Theta(n)$

并行算法计算复杂度 $Θ(n/p)\Theta(n/p)$ ，通信复杂度 $Θ(logp)\Theta(logp)$

并行开销 $To(n,p)=Θ(plogp)T_o(n,p) = \Theta(plogp)$

等效关系： $n≥Cplogpn\geq Cplogp$

问题规模为n时空间复杂度 $M(n)=Θ(n)M(n)=\Theta(n)$ ，即 $M (n) = n$

可扩展函数 $M (C p l o g p) / p = C p l o g p / p = C l o g p$

系统具有良好的可扩展性

第五章——MPI程序设计

在这里插入图片描述

MPI主要用于分布式内存并行系统。在分布式内存系统中，每个处理器/节点有自己独立的内存，不能直接访问其他处理器的内存。因此，不同进程之间需要通过消息传递来交换数据。MPI 正是为这种场景设计的一套标准通信接口，它允许开发者通过发送和接收消息在进程间进行通信，从而实现并行计算。

MPI并行程序中进程数量在启动时被指定且在整个程序执行过程中保持不变，每个进程都有唯一的ID号。MPI采用SIMD，所有进程执行相同的程序。

MPI函数

环境管理函数

1. MPI_Init

int MPI_Init(int *argc, char ***argv);

参数：

argc：指向 main 函数中 argc 的指针，用于传递命令行参数的数量。
argv：指向 main 函数中 argv 的指针，用于传递命令行参数的字符串数组。

返回值：成功时返回 MPI_SUCCESS，否则返回错误码。

MPI_Init 是 MPI 程序的第一个调用函数，用于初始化 MPI 库和并行执行环境，初始化后，所有进程默认属于 MPI_COMM_WORLD 通信域，这是一个包含所有进程的全局通信域，MPI_Init 必须与 MPI_Finalize 配对使用，后者用于释放 MPI 资源。

2. MPI_Finalize

int MPI_Finalize(void)

告知MPI系统MPI已经所使用完毕，为MPI而分配的所有资源都可以释放。

3. MPI_Comm_size

int MPI_Comm_size(MPI_Comm comm, int *size);

参数

comm：通信域（如 MPI_COMM_WORLD，表示所有进程的全局通信域）
size：输出参数，存储通信子中的进程数量

返回值：成功时返回 MPI_SUCCESS，否则返回错误码

返回指定通信域中的进程总数

4. MPI_Comm_rank

int MPI_Comm_rank(MPI_Comm comm, int *rank);

参数

comm：通信域
rank：输出参数，存储当前进程的 rank 值

返回值：成功时返回 MPI_SUCCESS，否则返回错误码

返回当前进程在通信域中的唯一标识符

集体通信

5. MPI_Reduce

int MPI_Reduce(
    void* sendbuf,         // 发送缓冲区（输入数据）
    void* recvbuf,         // 接收缓冲区（仅根进程有效）
    int count,             // 数据元素数量
    MPI_Datatype datatype, // 数据类型（如 MPI_INT, MPI_DOUBLE）
    MPI_Op op,             // 归约操作（如 MPI_SUM, MPI_MAX）
    int root,              // 结果接收进程的 rank
    MPI_Comm comm          // 通信域（如 MPI_COMM_WORLD）
);

所有进程的 sendbuf 数据通过 op 操作（如求和、求最大值）汇总到根进程的 recvbuf。采用树形网络通信，即第三章案例研究中规约问题的通信方式，通信复杂度为 $O (l o g P)$ 。

6. MPI_Allreduce

int MPI_Allreduce(
    void* sendbuf,         // 发送缓冲区（输入数据）
    void* recvbuf,         // 接收缓冲区（所有进程接收结果）
    int count,             // 数据元素数量
    MPI_Datatype datatype, // 数据类型（如 MPI_INT, MPI_DOUBLE）
    MPI_Op op,             // 归约操作（如 MPI_SUM, MPI_MAX）
    MPI_Comm comm          // 通信域（如 MPI_COMM_WORLD）
);

所有进程的 sendbuf 数据通过 op 操作（如求和、求最大值）聚合为一个结果，并存储于每个进程的 recvbuf。采用蝶形网络通信，蝶形网络每一轮通信次数多，但是发送方和接收方并不重叠，因此通信复杂度为 $O (l o g P)$ 。

在这里插入图片描述

7. MPI_Gather

int MPI_Gather(
    void* sendbuf,          // 发送缓冲区地址（每个进程需提供）
    int sendcount,          // 每个进程发送的数据量
    MPI_Datatype sendtype,  // 发送数据类型（如 MPI_INT）
    void* recvbuf,          // 接收缓冲区地址（仅根进程有效）
    int recvcount,          // 每个进程发送的数据量（需等于 sendcount）
    MPI_Datatype recvtype,  // 接收数据类型
    int root,               // 根进程的 rank
    MPI_Comm comm           // 通信域（如 MPI_COMM_WORLD）
);

每个进程（包括根进程）将自身 sendbuf 中的数据发送到根进程，根进程将这些数据按进程编号顺序拼接存储在 recvbuf 中。标准MPI库优先使用树形或蝶形算法，复杂度 $O (l o g P)$ ，仅在特定条件下退化为线性实现，复杂度 $O (P)$ 。

data = rank + 1;
int gathered_data[4];  // 根进程的接收缓冲区
MPI_Gather(&data, 1, MPI_INT, gathered_data, 1, MPI_INT, 0, MPI_COMM_WORLD);

if (rank == 0) {
    printf("Root collected: ");
    for (int i = 0; i < 4; i++) printf("%d ", gathered_data[i]);
    printf("\n");
}

Root collected: 1 2 3 4

变体函数

MPI_Gatherv：支持非均匀收集，允许不同进程发送不同数量的数据（需指定 displs 和 recvcounts 数组）。
MPI_Igather：非阻塞版本，支持异步通信。

8. MPI_Allgather

int MPI_Allgather(
    void* sendbuf,          // 发送缓冲区地址（每个进程提供）
    int sendcount,          // 每个进程发送的数据量
    MPI_Datatype sendtype,  // 发送数据类型（如 MPI_INT）
    void* recvbuf,          // 接收缓冲区地址（需预分配足够空间）
    int recvcount,          // 每个进程接收的数据量（通常等于 sendcount）
    MPI_Datatype recvtype,  // 接收数据类型
    MPI_Comm comm           // 通信域（如 MPI_COMM_WORLD）
);

将所有进程的数据收集并分发到通信域内的每个进程，使得每个进程最终都拥有全局数据的完整副本。主流MPI库采用树形结构优化，通信复杂度为 $O (l o g P)$ ，仅在特定条件下退化为线性实现，复杂度 $O (P)$ 。

变体函数

MPI_Allgatherv：支持非均匀数据收集（不同进程发送不同数据量）。
MPI_Iallgather：非阻塞版本，允许异步通信。

9. MPI_Scatter

int MPI_Scatter(
    void* sendbuf,          // 根进程的发送缓冲区（输入数据）
    int sendcount,          // 每个进程接收的数据量
    MPI_Datatype sendtype,  // 发送数据类型（如 MPI_INT）
    void* recvbuf,          // 接收缓冲区（存储分发到的数据）
    int recvcount,          // 接收的数据量（通常等于 sendcount）
    MPI_Datatype recvtype,  // 接收数据类型
    int root,               // 根进程的 rank
    MPI_Comm comm           // 通信域（如 MPI_COMM_WORLD）
);

根进程（root）将一个数组的不同部分分发给各进程，按进程编号顺序分配。MPI_Scatter的标准实现是树形分发，通信复杂度为 $O (l o g P)$ ，若是线性分发，通信复杂度为 $O (P)$ 。MPI_Scatter适合块划分场景。

int send_data[6] = {1, 2, 3, 4, 5, 6};
int recv_data[2];  
MPI_Scatter(send_data, 2, MPI_INT, recv_data, 2, MPI_INT, 0, MPI_COMM_WORLD);
printf("Process %d received: %d, %d\n", rank, recv_data[0], recv_data[1]);

Process 0 received: 1, 2
Process 1 received: 3, 4
Process 2 received: 5, 6

变体函数

MPI_Scatterv：支持非均匀分发，允许不同进程接收不同数量的数据（需指定偏移量数组 displs）。
MPI_Iscatter：非阻塞版本，允许异步通信。

10. MPI_Bcast

int MPI_Bcast(
    void* buffer,          // 发送/接收缓冲区指针（根进程发送，其他进程接收）
    int count,             // 数据元素个数
    MPI_Datatype datatype, // 数据类型（如 MPI_INT、MPI_DOUBLE）
    int root,              // 根进程的 rank
    MPI_Comm comm          // 通信域（如 MPI_COMM_WORLD）
);

根进程的缓冲区数据会被复制到所有其他进程的缓冲区中，确保所有进程获得相同数据。主流MPI库使用树形结构优化广播效率，通信复杂度 $O (l o g P)$ 。MPI_Bcast适合所有进程需要相同数据的初始化数据分发。

同步控制函数

11. MPI_Barrier

int MPI_Barrier(MPI_Comm communicator);

阻塞调用进程，直到通信域内所有进程都调用 MPI_Barrier，形成逻辑上的“屏障”。由于集体通信隐含了同步机制，因此集体通信无需显式调用 MPI_Barrier 。

应用场景

计时同步

MPI_Barrier(MPI_COMM_WORLD);
double start = MPI_Wtime();
// 并行计算代码
MPI_Barrier(MPI_COMM_WORLD);
double end = MPI_Wtime();

阶段控制：在分布式计算中划分阶段，例如确保所有进程完成数据加载后再开始计算
调试辅助：强制进程同步以定位执行不一致的问题

计时函数

12. MPI_Wtime

double MPI_Wtime(void);

返回一个双精度浮点数，表示从某个固定时刻（如程序启动或MPI初始化）到调用时刻的时间。比标准C的clock()更准确，后者仅测量CPU时间（不包括等待时间）。

13. MPI_Wtick

double MPI_Wtick(void);

返回MPI_Wtime的时间分辨率（即时钟滴答周期），单位为秒。

double tick = MPI_Wtick();
printf("时间精度: %e 秒\n", tick);

点对点通信

14. MPI_Send

int MPI_Send(
    void* buf, 				// 发送数据的缓冲区地址
    int count, 				// 发送的数据项数量
    MPI_Datatype datatype,	// 数据类型（如MPI_INT）
    int dest, 				// 目标进程的rank（编号）
    int tag, 				// 消息标签，用于区分不同消息
    MPI_Comm comm			// 通信域（如MPI_COMM_WORLD）
);

该函数会阻塞直到直到消息缓冲区空闲，消息缓冲区空闲的时间是消息被复制到系统缓冲区或传送，典型情况下消息被复制到系统缓冲区，传输重叠的计算。

15. MPI_Recv

int MPI_Recv(
    void* buf, 				// 发送数据的缓冲区地址
    int count, 				// 发送的数据项数量
    MPI_Datatype datatype,	// 数据类型（如MPI_INT）
    int source, 			// 发送进程的rank，可用MPI_ANY_SOURCE接收任意进程的消息。
    int tag, 				// 消息标签，用于区分不同消息
    MPI_Comm comm, 			// 通信域（如MPI_COMM_WORLD）
    MPI_Status* status		// 返回接收状态（如消息来源和实际接收的数据量）
);

函数阻塞，直到信息进入缓冲区，如果信息没有到达，函数就不会返回。

死锁

死锁:等待一个永远不会成真的条件的过程

容易写出死锁的发送/接收代码：

两个进程：都是先收后发
发送标签与接收标签不一致
进程向错误的目标进程发送消息

16. MPI_Ssend

int MPI_Ssend(
    void* buf,          	// 发送缓冲区地址
    int count,         		// 数据项数量
    MPI_Datatype datatype, 	// 数据类型（如 MPI_INT）
    int dest,           	// 目标进程的 rank
    int tag,            	// 消息标签（用于区分消息）
    MPI_Comm comm       	// 通信域（如 MPI_COMM_WORLD）
);

特性	`MPI_Send`	`MPI_Ssend`
同步要求	可能缓存数据或直接发送	必须等待接收方启动接收
返回条件	数据被缓存或发送完成	接收方已开始接收且数据已发送完成
性能	通常更快（允许系统优化）	较慢（需同步握手）
安全性	可能因未匹配接收导致死锁	更易避免死锁（强制同步）

17. MPI_Sendrecv

int MPI_Sendrecv(
    void* sendbuf,         // 发送缓冲区地址
    int sendcount,         // 发送数据项数量
    MPI_Datatype sendtype, // 发送数据类型（如 MPI_INT）
    int dest,              // 目标进程 rank
    int sendtag,           // 发送消息标签
    void* recvbuf,         // 接收缓冲区地址
    int recvcount,         // 接收数据项数量
    MPI_Datatype recvtype, // 接收数据类型
    int source,            // 源进程 rank（可设为 MPI_ANY_SOURCE）
    int recvtag,           // 接收消息标签（可设为 MPI_ANY_TAG）
    MPI_Comm comm,         // 通信域（如 MPI_COMM_WORLD）
    MPI_Status* status     // 返回接收状态信息
);

阻塞式操作：函数会阻塞，直到发送完成且接收完成后才返回，发送和接收使用同一通信域，但可指定不同的标签（sendtag 和 recvtag）。
缓冲区分离：发送缓冲区（sendbuf）和接收缓冲区（recvbuf）必须独立，且可具有不同长度和数据类型；
避免死锁：在环形通信或链式通信中（如进程间两两交换数据），MPI_Sendrecv 能自动处理依赖关系，无需手动调整发送/接收顺序；
兼容性：可与普通MPI_Send、MPI_Recv混合使用，即发送的消息可被普通接收操作接收，反之亦然。

分配方式

按我的理解，分配方式体现了Foster设计方法的聚集和映射阶段。

交错分配/循环分配

交替分配/循环分配均通过交替或循环的方式将任务分配给不同进程。

循环分配：每个进程按固定步长（如进程数p）循环获取任务。
交错分配：任务按交替顺序分配给进程（如轮询分配），本质与循环分配类似。

如电路可满足问题中

for (i = id; i < 65536; i += p) {
    count += check_circuit(id, i);
}

任务 $i$ 通过模运算 $pi\mod p$ 聚集到不同的逻辑组 $i d$ 中，然后将聚集后逻辑组 $i d$ 的映射到具体的物理进程 $i d$ 。

块分配

块分配是一种常见的数据划分策略，用于将任务或数据均匀分配给多个处理单元。

如在素数筛选中

low_value = 2 + id * (n - 1) / p;
high_value = 1 + (id + 1) * (n - 1) / p;
size = high_value - low_value + 1;
marked = (char*)malloc(size);

将数字范围 [2, n] 连续划分为 p 个块（p 为进程数），每个进程 id 负责处理一个连续子区间 [low_value, high_value]。

[2, n] 的最小划分是 n-1 个数字，这里将其聚集为 p 个相对均匀的块，因为是根据每个进程的 $i d$ 计算出，所以相当于将这 p 各块映射到 p 个处理单元。

电路可满足性问题

在这里插入图片描述
电路可满足性问题是NP完全的，没有已知的算法可以在多项式时间内解决，因此我们寻找所有的解决方案，通过穷举式搜索16个输入共65,536个组合来测试。

并行算法设计

划分：可能的65,536组合，检查每个组合是基本任务
通信：基本任务之间无需通信；各进程统计本地合法组合数量后，通过 MPI_Reduce 汇总到进程0
聚集：按组合可能 $i$ ， $pi\mod p$ 逻辑分组
映射：绑定逻辑分组到进程 $i d$

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

#define EXTRACT_BIT(n, i) ((n & (1 << i)) ? 1 : 0)

int check_circuit (int id, int z) {
    int v[16];
    int i;
    
    for (i = 0; i < 16; i++) v[i] = EXTRACT_BIT(z, i);
    
    if ((v[0] || v[1]) && (!v[1] || !v[3]) && (v[2] || v[3])
        && (!v[3] || !v[4]) && (v[4] || ！v[5])
        && (v[5] || !v[6]) && (v[5] || v[6])
        && (v[6] || !v[15]) && (v[7] || !v[8])
        && (!v[7] || !v[13]) && (v[8] || v[9])
        && (v[8] || !v[9]) && (!v[9] || !v[10])
        && (v[9] || v[11]) && (v[10] || v[11])
        && (v[12] || v[13]) && (v[13] || !v[14])
        && (v[14] || v[15])) {
        printf("%d) %d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d\n", id,
              v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7],
              v[8], v[9], v[10], v[11], v[12], v[13], v[14], 
              v[15],
        );
        fflush(stdout);
        return 1;
    }
    return 0;
}

int main (int argc, char *argv[]) {
    int i;
    int id;
    int p;
    int count;
    int global_count;
    double elapsed_time;
    
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    MPI_Barrier (MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();
    
    count = 0;
    for (i = id; i < 65536; i += p) {
        count += check_circuit(id, i);
    }
    
    MPI_Reduce(&count, &global_count, 1, MPI_INT,
               MPI_SUM, 0, MPI_COMM_WORLD);
    
	MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time += MPI_Wtime();
    
    printf("Process %d is done, cost %lf sec\n", id, elapsed_time);
    fflush(stdout);
    
    if (!id) 
        printf("There are %d different solutions\n", global_count);
        
    MPI_Finalize();
    return 0;
}

素数筛选——埃拉托斯特尼筛法

在这里插入图片描述

串行算法

创建无标记的自然数列表 2, 3, …, n
$k = 2$
重复进行
- (a) 标记 $k^{2}$ 到 $n$ 之间 $k$ 的所有倍数
- (b) 将最小的未标记的且大于 $k$ 的数赋值给 $k$ *，直到 $k^{2} > n$
未标记的数就是素数

外层循环的 $p$ 从 2 到 $n\sqrt{n}$ ，因此循环次数约为 $O(n)O(\sqrt{n})$ ；每个素数 $p$ 会标记 $p^2, p^2 + p, p^2 + 2p, ..., n$ 中的所有倍数，对于每个素数 $p$ ，标记次数大约是 $⌊n−p2p⌋+1≈np\left\lfloor \frac{n - p^2}{p} \right\rfloor + 1 \approx \frac{n}{p}$ ；所以总标记操作的次数近似为：
$\sum_{p \text{为素数}} \frac{n}{p} \approx n \cdot \sum_{p \le \sqrt{n}} \frac{1}{p}$
由于调和级数中素数倒数和的渐近性质为 $∑p≤x1p∼ln⁡ln⁡x\sum_{p \le x} \frac{1}{p} \sim \ln \ln x$ ，我们有：
$\text{总标记操作数} \approx n \cdot \ln \ln n$
假设每次标记耗时 $μ\mu$ ，总耗费时间为 $T(n)=μ∗n∗ln(ln(n))T(n)=\mu * n * ln(ln(n))$ 。

并行算法设计

划分：将筛选范围 [2, n] 的数字划分为独立的标记任务，每个数字的标记操作是天然并行的，只需知道当前素数 prime 即可标记其倍数
通信：进程0负责找到下一个素数 prime，通过 MPI_Bcast 广播给所有进程；各进程统计本地素数数量后，通过 MPI_Reduce 汇总到进程0
聚集：将连续的数字区间 [low_value, high_value] 合并为粗粒度任务块
映射：进程 id 固定负责区间 [low_value, high_value] 的标记操作

#include "mpi.h"
#include <math.h>
#include <stdio.h>
#include<stdlib.h>

#define MIN(a,b)  ((a)<(b)?(a):(b))
#define LL long long

int main(int argc, char* argv[]){
    LL    	count;        		/* Local prime count */
    double 	elapsed_time; 		/* Parallel execution time */
    LL    	first;        		/* Index of first multiple */
    LL    	global_count = 0; 	/* Global prime count */
    LL    	high_value;   		/* Highest value on this proc */
    LL    	i;
    int    	id;           		/* Process ID number */
    LL    	index;        		/* Index of current prime */
    LL    	low_value;    		/* Lowest value on this proc */
    char* 	marked;       		/* Portion of 2,...,'n' */
    LL    	n;            		/* Sieving from 2, ..., 'n' */
    int    	p;           		/* Number of processes */
    LL    	proc0_size;   		/* Size of proc 0's subarray */
    LL    	prime;        		/* Current prime */
    LL    	size;         		/* Elements in 'marked' */

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);
    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time = -MPI_Wtime();

    if (argc != 2) {
        if (!id) printf("Command line: %s <m>\n", argv[0]);
        MPI_Finalize();
        exit(1);
    }

    n = atoll(argv[1]);

    proc0_size = (n - 1) / p;
    if ((2 + proc0_size) < (int)sqrt((double)n)) {
        if (!id) printf("Too many processes\n");
        MPI_Finalize();
        exit(1);
    }

    low_value = 2 + id * (n - 1) / p;
    high_value = 1 + (id + 1) * (n - 1) / p;
    size = high_value - low_value + 1;

    marked = (char*)malloc(size);
    if (marked == NULL) {
        printf("Cannot allocate enough memory\n");
        MPI_Finalize();
        exit(1);
    }
    
    for (i = 0; i < size; i++) marked[i] = 0;
    
    if (!id) index = 0;
    prime = 2;
    do {
        if (prime * prime > low_value)
            first = prime * prime - low_value;
        else {
            if (!(low_value % prime)) first = 0;
            else first = prime - (low_value % prime);
        }
        for (i = first; i < size; i += prime) marked[i] = 1;
        if (!id) {
            while (marked[++index]);
            prime = index + 2;
        }
        if (p > 1) MPI_Bcast(&prime, 1, MPI_INT, 0, MPI_COMM_WORLD);
    } while (prime * prime <= n);

    count = 0;
    for (i = 0; i < size; i++) if (!marked[i]) count++;
    MPI_Reduce(&count, &global_count, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);

    MPI_Barrier(MPI_COMM_WORLD);
    elapsed_time += MPI_Wtime();

    if (!id) {
        printf("There are %lld primes less than or equal to %lld\n", global_count, n);
        printf("SIEVE (%d) %10.6f\n", p, elapsed_time);
    }
    MPI_Finalize();
    return 0;
}

分析

标记时间由 $p$ 个进程共享，耗时为 $T1(n)=μ∗np∗ln(ln(n))T_1(n)=\mu * \frac{n}{p} * ln(ln(n))$ ；广播一个整数，耗时近似为 $λ\lambda$ ，一共要广播大约 $π(n)≈nln⁡n\pi(\sqrt{n}) \approx \frac{\sqrt{n}}{\ln \sqrt{n}}$ 个素数，每次广播需要 $⌈logp⌉\lceil logp \rceil$ 次通信，因此耗时为 $T2=λ∗nln⁡n∗⌈logp⌉T_2= \lambda * \frac{\sqrt{n}}{\ln \sqrt{n}} * \lceil logp \rceil$ 。总耗时：
$T_n = \mu * \frac{n}{p} * ln(ln(n)) + \lambda * \frac{\sqrt{n}}{\ln \sqrt{n}} * \lceil logp \rceil$
优化方法

删除偶数：将计算次数减半，为更大的n值释放存储空间
每个进程找到自己的筛选素数：将素数的计算复制到 $n\sqrt n$ ，消除了广播步骤
重新组织循环：增加缓存命中率

弗洛伊德算法

串行算法

Input: 
    n             // 顶点个数
    D[1..n][1..n] // 邻接矩阵，D[i][j] 表示从顶点 i 到 j 的初始距离
                  // 若 i == j 则 D[i][j] = 0，若 i ≠ j 且无边则 D[i][j] = ∞

Algorithm:
    for k from 1 to n do
        for i from 1 to n do
            for j from 1 to n do
                if D[i][j] > D[i][k] + D[k][j] then
                    D[i][j] := D[i][k] + D[k][j]

Output: 
    D[1..n][1..n] // 所有点对的最短路径长度

显然时间复杂度为 $O(n^{3})$ ，若判断加标记时间为 $μ\mu$ ，则时耗为 $\mu * n^{3}$ 。

并行算法设计

划分：串行算法中并没有功能并行，因此选择数据划分，列块方式划分消除了列的广播，行块方式划分消除了行内广播并且从文件中读取矩阵更简单，因此选择行块方式划分。
通信：Floyd 算法中每行是独立更新的单位，D[i][j] 的更新只依赖于 D[i][j] 和 D[i][k], D[k][j]，对于固定的 i，整行 D[i][*] 是一起更新的。如果一个进程拥有第 i 行的数据，就可以本地完成这行的全部计算（只需要 k 行作为中间数据）。因此每一轮 k，需要将“第 k 行”广播给所有进程。
聚集：行块划分 + 同步广播，避免了每次都单独通信每个元素（用广播一次性发送整行），将较小任务“聚合”为一批元素处理（一个进程一次处理多行）。
映射：使用了 BLOCK_SIZE、BLOCK_LOW、BLOCK_OWNER 等宏，确保每个进程获得相近的任务量，BLOCK_OWNER(k, p, n) 负责确定每行属于哪个进程，实现负载均衡的静态映射。

分析

因为将行进行块划分，因此本地计算耗时 $T1=μ∗n∗⌈np⌉∗n=μ∗n2∗⌈np⌉T_1=\mu * n * \lceil \frac{n}{p} \rceil *n = \mu * n^{2} * \lceil \frac{n}{p} \rceil$ ，总共需要进行k次广播，每次广播共有 $⌈logp⌉\lceil logp \rceil$ 次通信，每次通信耗时 $λ+4nβ\lambda + \frac{4n}{\beta}$ ，其中 $λ\lambda$ 是传播时间， $4nβ\frac{4n}{\beta}$ 是传输时间，通信耗时 $T2=n⌈logp⌉(λ+4nβ)T_2=n\lceil logp \rceil(\lambda + \frac{4n}{\beta})$ 。总耗时：
$\mu * n^{2} * \lceil \frac{n}{p} \rceil + n\lceil logp \rceil(\lambda + \frac{4n}{\beta})$

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <limits.h>

#define dtype int
#define MPI_TYPE MPI_INT
#define MIN(a, b) ((a) < (b) ? (a) : (b))

#define BLOCK_LOW(id, p, n) ((id)*(n)/(p))
#define BLOCK_HIGH(id, p, n) (BLOCK_LOW((id)+1, p, n) - 1)
#define BLOCK_SIZE(id, p, n) (BLOCK_HIGH(id, p, n) - BLOCK_LOW(id, p, n) + 1)
#define BLOCK_OWNER(k, p, n) (((p)*(k)+n-1)/(n))

void terminate(int id, const char *message) {
    if (id == 0) fprintf(stderr, "%s", message);
    MPI_Abort(MPI_COMM_WORLD, 1);
}

void read_row_striped(char *filename, void ***a, void **storage, 
                      MPI_Datatype dtype, 
                      int *m, int *n, 
                      MPI_Comm comm) {
    int id, p;
    FILE *infile;
    int i, j;
    int local_rows;

    MPI_Comm_rank(comm, &id);
    MPI_Comm_size(comm, &p);

    if (id == 0) {
        infile = fopen(filename, "r");
        if (!infile) terminate(id, "Cannot open file\n");
        fscanf(infile, "%d %d", m, n);
    }

    MPI_Bcast(m, 1, MPI_INT, 0, comm);
    MPI_Bcast(n, 1, MPI_INT, 0, comm);

    local_rows = BLOCK_SIZE(id, p, *m);
    *storage = malloc(local_rows * (*n) * sizeof(dtype));
    *a = malloc(local_rows * sizeof(dtype *));

    for (i = 0; i < local_rows; i++)
        (*a)[i] = &((dtype *)(*storage))[i * (*n)];

    if (id == 0) {
        dtype *full_mat = malloc((*m) * (*n) * sizeof(dtype));
        for (i = 0; i < *m; i++)
            for (j = 0; j < *n; j++)
                fscanf(infile, "%d", &full_mat[i * (*n) + j]);
        fclose(infile);

        int dest;
        for (dest = 0; dest < p; dest++) {
            int low = BLOCK_LOW(dest, p, *m);
            int rows = BLOCK_SIZE(dest, p, *m);
            if (dest == 0) {
                for (i = 0; i < rows; i++)
                    for (j = 0; j < *n; j++)
                        (*a)[i][j] = full_mat[(low + i) * (*n) + j];
            } else {
                MPI_Send(&full_mat[low * (*n)], rows * (*n), dtype, dest, 0, comm);
            }
        }
        free(full_mat);
    } else {
        MPI_Recv(*storage, local_rows * (*n), dtype, 0, 0, comm, MPI_STATUS_IGNORE);
    }
}

void print_row_striped(void **a, 
                       MPI_Datatype dtype, 
                       int m, int n, 
                       MPI_Comm comm) {
    int id, p;
    MPI_Comm_rank(comm, &id);
    MPI_Comm_size(comm, &p);

    int i, j;
    void *buf = NULL;

    if (id == 0) {
        buf = malloc(n * sizeof(dtype));
        for (i = 0; i < m; i++) {
            int owner = BLOCK_OWNER(i, p, m);
            if (owner == 0) {
                for (j = 0; j < n; j++)
                    printf("%d ", ((dtype **)a)[i][j]);
                printf("\n");
            } else {
                MPI_Recv(buf, n, dtype, owner, 0, comm, MPI_STATUS_IGNORE);
                for (j = 0; j < n; j++)
                    printf("%d ", ((dtype *)buf)[j]);
                printf("\n");
            }
        }
        free(buf);
    } else {
        int local_low = BLOCK_LOW(id, p, m);
        int local_high = BLOCK_HIGH(id, p, m);
        int rows = local_high - local_low + 1;
        for (i = 0; i < rows; i++) {
            MPI_Send(((dtype **)a)[i], n, dtype, 0, 0, comm);
        }
    }
}

void compute_shortest_paths(int id, int p, dtype **a, int n) {
    int i, j, k;
    int offset;
    int root;
    dtype *tmp;

    tmp = (dtype *)malloc(n * sizeof(dtype));
    for (k = 0; k < n; k++) {
        root = BLOCK_OWNER(k, p, n);
        if (root == id) {
            offset = k - BLOCK_LOW(id, p, n);
            for (j = 0; j < n; j++)
                tmp[j] = a[offset][j];
        }

        MPI_Bcast(tmp, n, MPI_TYPE, root, MPI_COMM_WORLD);

        for (i = 0; i < BLOCK_SIZE(id, p, n); i++) {
            for (j = 0; j < n; j++) {
                a[i][j] = MIN(a[i][j], a[i][k] + tmp[j]);
            }
        }
    }
    free(tmp);
}

int main(int argc, char *argv[]) {
    dtype **a;
    dtype *storage;
    int i, j, k;
    int id;
    int m;
    int n;
    int p;
    double time, max_time;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &p);

    read_row_striped(argv[1], (void *)&a, (void *)&storage, 
                     MPI_TYPE, &m, &n, MPI_COMM_WORLD);
    if (m != n) terminate(id, "Matrix must be square\n");

    print_row_striped((void **)a, MPI_TYPE, m, n, MPI_COMM_WORLD);

    MPI_Barrier(MPI_COMM_WORLD);
    time = -MPI_Wtime();

    compute_shortest_paths(id, p, (dtype **)a, n);

    time += MPI_Wtime();
    MPI_Reduce(&time, &max_time, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

    if (id == 0)
        printf("Floyd, matrix size %d, %d processes: %.6f seconds\n", n, p, max_time);
        
    print_row_striped((void **)a, MPI_TYPE, m, n, MPI_COMM_WORLD);
    free(a);
    free(storage);
    MPI_Finalize();
    return 0;
}

标签

#笔记 #算法 #c语言