I.3 科学应用的特征

可扩展共享内存多处理器的主要用途是真正的并行程序设计，而非多道程序运行或面向事务的计算。并行计算主要面向科学与工程技术应用。因此，要理解其设计问题，需要先了解这类应用的行为。本节对此作一个介绍。

科学应用的特征

本节的科学与工程并行工作负载由两个应用和两个计算内核组成。两个内核分别是快速傅里叶变换（FFT）和 LU 分解。之所以选择它们，是因为它们代表了广泛应用中常用的技术，而且其性能特征具有许多并行科学应用的典型特点。此外，这两个内核的代码较短，容易理解其行为，并直接追溯到具体的体系结构特征。和许多科学应用一样，这组工作负载基本没有 I/O。

本附录使用的两个应用是 Barnes 和 Ocean，代表了两种重要、但很不相同的并行计算类型。下面简要介绍各应用和内核，从并行度与通信两个方面刻画其基本行为。我们将说明如何为分布式共享内存多处理器分解问题；对于只使用一套集中式存储器的多处理器，下面某些数据分解并非必要。

FFT 内核

FFT 是采用谱方法的应用中的关键内核，涉及的领域从信号处理、流体流动到气候建模。本节研究的 FFT 应用，是复数 FFT 并行算法的一维版本。对 $n$ 个数据点，其顺序执行时间为 $n\log n$。算法采用等于 $\sqrt n$ 的高基数，以尽量减少通信。本附录的测量采用含一百万个点的输入数据集。

算法有三个主要数据结构：待变换数据的输入数组、输出数组，以及预先计算、执行期间只读的单位根矩阵。所有数组均组织成方阵。算法有六个步骤：

1. 转置数据矩阵。
2. 对数据矩阵的每一行执行一维 FFT。
3. 将单位根矩阵与数据矩阵相乘，结果写回数据矩阵。
4. 转置数据矩阵。
5. 对数据矩阵的每一行执行一维 FFT。
6. 转置数据矩阵。

数据矩阵和单位根矩阵按连续的若干行划分给各处理器，使各处理器负责的部分位于其本地存储器。所有处理器都会频繁访问单位根矩阵的第一行，因此常在上述第一步复制这一行；这里也采用这种做法。数据转置保证了各 FFT 步骤有良好的局部性，否则这些步骤会访问非本地数据。

通信只发生在转置阶段，需要进行大量的全到全通信（all-to-all）。将分配给一个处理器的若干行中的连续列片段组成块，转置这些块，并将其放到目的矩阵中的正确位置。每个处理器在本地转置一个块，同时向系统中其他每个处理器各发送一个块。虽然转置中没有对单个字的重用，但缓存块较长时，仍值得将转置分块，以利用源矩阵中长缓存块所提供的空间局部性。

LU 内核

LU 对稠密矩阵作 LU 分解，代表了 QR 分解、Cholesky 分解和特征值方法等许多稠密线性代数计算。对一个 $n\times n$ 矩阵，运行时间为 $n^3$，并行度与 $n^2$ 成正比。利用第 2 章介绍的技术对算法分块，可以高效执行稠密 LU 分解，获得很好的缓存行为，并减少通信。分块之后，主要计算是最内层循环中的稠密矩阵乘法。块应足够小，以保持较低的缓存未命中率；同时也应足够大，以减少在计算中并行度较低部分所花的时间。较小的块，例如 $8\times8$ 或 $16\times16$，通常就能同时满足这两个条件。

减少处理器间通信有两个重要细节。首先，以二维平铺方式将矩阵块分配给处理器：对于由 $\frac nB\times\frac nB$ 个块组成的矩阵，每块为 $B\times B$，将一个 $p\times p$ 网格反复铺在块矩阵上，直到全部块都分配给处理器。其次，稠密矩阵乘法由拥有目的块的处理器完成。采用这种分块与分配方案后，消元过程中的通信既规则又可预测。本附录的测量使用 $512\times512$ 输入矩阵和 $16\times16$ 的块。

在共享地址空间中编写二维矩阵的分块 LU 分解，一种自然方法是用二维数组表示矩阵。但由于块按平铺方式分配，而一个块在二维数组的地址空间中并不连续，因此很难把各块安置在拥有它们的处理器的本地存储器中。解决办法是使用四维数组，保证分配给同一处理器的块连续地存放在本地：前两维指定块在二维块网格中的编号，后两维指定块内元素。

译注： 上段网格大小为 $p\times p$，其边长符号不能直接当成图 I.1 中表示处理器总数的 $p$，否则会把网格位置数与处理器数量混为一谈。

Barnes 应用

Barnes 实现了 Barnes–Hut n 体算法，用于求解星系演化中的一个问题。n 体算法模拟大量物体之间的相互作用力。在这里，物体代表恒星的集合，相互作用力是引力。完整计算全部物体两两之间的相互作用，计算量按 $n^2$ 增长。为降低计算时间，n 体算法利用了力随距离衰减这一性质。例如，引力按 $1/d^2$ 衰减，其中 $d$ 为两物体间的距离。Barnes–Hut 算法把距离某个物体“足够远”的一组物体，近似为位于该组质心处的一个点，其质量等于整组物体的质量。如果目标物体与该组中任一物体都足够远，这种近似引入的误差就可以忽略。各集合按层次组织，可以用树表示。算法的运行时间为 $n\log n$，并行度与 $n$ 成正比。

Barnes–Hut 算法使用八叉树，每个节点最多有八个子节点，用来表示一部分空间中的八个立方体。每个节点代表以其为根的子树中的全部物体，我们把这一集合称为空间单元（cell）。由于空间密度不同，而叶节点代表单个物体，树的深度也不同。对每个物体，都遍历一次树，计算作用于该物体的合力。力的计算从树根开始：每访问一个节点，判断以该节点为根的子树所表示空间单元的质心，离当前物体是否“足够远”。如果足够远，就以位于该单元质心的一个点近似整棵子树，计算这个质心对该物体施加的力；否则，必须“打开”该单元，访问它的各棵子树。物体与单元的距离，以及允许的误差，共同决定哪些单元需要打开。计算力的阶段占据主要执行时间。本附录用 16K 个物体进行测量，判断单元是否需要打开的参数取实际常用范围的中间值。

Barnes–Hut 要获得高效并行性能并不容易。物体分布不均匀，而且随时间变化，因而难以在处理器间划分工作，并保持良好的引用局部性。两个性质对此有帮助：第一，系统演化缓慢；第二，引力随距离迅速衰减，所以每个单元很可能只需访问少数其他单元，其中大部分在上一时间步已经访问过。可以通过给每个处理器分配一棵子树来划分整棵树。计算子树中一个物体所受力时，许多访问仍然指向同一子树中的其他物体。由于各子树的工作量不同——空间中密集区域的单元需要访问更多单元——分配给处理器的子树大小，应依据某种工作量度量，例如需要访问多少其他单元，而不只是依据子树的节点数。对八叉树表示进行这种划分，可以兼顾负载平衡和引用局部性，同时保持较低的划分成本。不过，尽管局部性很好，最终的数据访问仍往往规模很小、结构不规则。因此，这种方案需要高效的共享内存通信实现。

Ocean 应用

Ocean 模拟涡流和边界流对大尺度海洋流动的影响。它用受限的红黑 Gauss–Seidel 多重网格技术，求解一组椭圆型偏微分方程。红黑 Gauss–Seidel 是一种迭代技术：给网格点着色，从而能够一致地依据相邻点的旧值更新各点。多重网格方法使用分层网格，通过迭代求解有限差分方程。层次中的每一层网格都比下一层点数更少，是对下一层网格的近似。较细的网格提高精度，从而提高收敛速率，但因数据点更多，也需要更多执行时间。下一次迭代应沿网格层次向上还是向下移动，取决于数据值的变化速率。每个时间步的误差估计用于决定继续留在当前网格、转向更粗网格，还是转向更细网格。当迭代在最细层次收敛时，就得到了一个解。对于 $n\times n$ 网格，每次迭代的工作量为 $n^2$，并行度也具有相同规模。

表示各层网格的数组根据具体问题动态分配并确定大小。整个海盆被划分为尽可能接近正方形的子网格，分配到各处理器本地存储器所对应的地址空间中；各处理器负责自己的子网格。本附录测量使用含 $130\times130$ 个网格点的输入。一次时间迭代分为五步。由于各步之间需要交换数据，所有处理器在每一步结束时同步，再进入下一步。当相邻子网格按最近邻方式访问彼此边界上的点时，就会发生通信。

并行程序的计算与通信

决定并行程序性能的一个关键特征，是计算量与通信量的比值。如果该比值很高，就意味着每传送一个数据项，应用会进行大量计算。如 I.2 节所述，通信是并行计算中代价高昂的部分，因此较高的计算／通信比十分有利。在并行处理环境中，我们关心的是：增加处理器数量、扩大问题规模，或同时改变二者时，这个比值如何变化。了解它随处理器数增长的变化，有助于判断应用能够获得多大加速。由于我们也常希望求解更大的问题，因此了解数据集大小对这一比值的影响同样关键。

要定量理解增加处理器后的计算／通信比，可以分别考察增加处理器或扩大问题规模时，计算量与通信量如何变化。图 I.1 表明，对这些应用，随着处理器数量增加，每个处理器的计算量按比例下降，而每个处理器的通信量下降得较慢。问题规模扩大时，计算量按算法的 $O(\cdot)$ 复杂度增长。通信的变化更复杂，取决于算法细节；图 I.1 的说明概括了各应用的基本规律。

计算量随 $n$ 按 $O(\cdot)$ 分析给出的速率增长，随 $p$ 增加则线性下降。通信的变化更复杂。在 FFT 中，所有数据点都必须相互作用，所以通信量随 $n$ 增加、随 $p$ 减少。在 LU 和 Ocean 中，通信量与块的边界成正比，因此随数据集大小增长的速率，与包含 $n$ 个点的正方形边长——即 $\sqrt n$——成正比；同理，二者的通信量与 $\sqrt p$ 成反比。Barnes 最为复杂：由于物体间相互作用随距离衰减，需要通信的基本相互作用数量按 $\sqrt n$ 增长；为维护物体之间的关系，还需要一个 $\log n$ 因子。处理器数增加时，通信量按 $1/\sqrt p$ 变化。

译注： 表中 LU 的计算量按原表保留为 $n/p$。若这里的 $n$ 按图注解释为矩阵元素总数，则完整稠密 LU 分解的计算量与前文的立方阶描述并不一致。原书没有在此交代额外归一化或单步计算假设，不据此改造原表。

整体计算／通信比由计算量与通信量各自的增长速率求得。通常，这个比值随数据集增大而缓慢上升，随处理器数增加而下降。这提醒我们：用更多处理器求解固定规模问题时，处理器间的通信增长，会带来越来越多的低效。它还说明，为使通信时间所占比例不变，增加处理器时，数据集需要以多快的速度增长。下面的例题展示这些取舍。

例题 假设已知 Ocean 在某台多处理器上使用四个处理器运行时，执行时间的 20% 用于等待通信。假定每次通信事件的代价与处理器数无关；一般情况下并非如此，因为通信代价会随处理器数增加。对同一规模的问题，预计在 32 处理器机器上能快多少？此时通信占执行时间的多少？若想保持通信时间占比不变，应将问题扩大多少？

原书解答 Ocean 的计算／通信比为 $\sqrt n/\sqrt p$。因此，问题规模不变时，通信频率按 $\sqrt p$ 增长。原书据此认为通信时间增加到 $\sqrt8$ 倍，并用 Amdahl 定律的一种变形，同时计入计算时间减少和通信时间增加。设四个处理器的总执行时间为 $T_4$，则 32 个处理器的执行时间为：

于是加速比为：

通信时间占比从 20% 变成 $0.57/0.67=85\%$。

为了使通信时间占比保持不变，计算／通信比必须保持不变，因此问题规模必须与处理器数以相同的速率增长。注意，此时问题规模已经改变，不能再公平地比较原问题与扩大后问题的加速比。I.6 节会再次讨论针对多处理器扩大应用规模这一关键问题。

译注： 上述解答与图 I.1 的“每处理器”通信量不一致。按表中 Ocean 的 $\sqrt n/\sqrt p$ 以及题设通信单价不变的假定，固定问题由 4 扩到 32 个处理器时，每处理器通信时间应除以 $\sqrt8$，不是乘以 $\sqrt8$。相应地，$T_{32}/T_4=0.8/8+0.2/\sqrt8\approx0.17071$，加速比约 5.86，通信占比约 41.4%。原书把相对于计算的通信比例上升，直接用成了绝对通信时间上升。这里完整保留原解答，并明确区分这一按题设计算的校正；保持通信占比需将数据集扩大 8 倍的结论不受此处方向错误影响。