# 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 的说明概括了各应用的基本规律。 | 应用 | 计算量的变化 | 通信量的变化 | 计算／通信比的变化 | | --- | --- | --- | --- | | FFT | \(\dfrac{n\log n}{p}\) | \(\dfrac np\) | \(\log n\) | | LU | \(\dfrac np\) | \(\dfrac{\sqrt n}{\sqrt p}\) | \(\dfrac{\sqrt n}{\sqrt p}\) | | Barnes | \(\dfrac{n\log n}{p}\) | 约为 \(\dfrac{\sqrt n\log n}{\sqrt p}\) | 约为 \(\dfrac{\sqrt n}{\sqrt p}\) | | Ocean | \(\dfrac np\) | \(\dfrac{\sqrt n}{\sqrt p}\) | \(\dfrac{\sqrt n}{\sqrt p}\) | {#fig-I.1 .wide .expand num="I.1" caption="计算量、通信量及二者比值的扩展规律，是决定并行多处理器性能的关键因素。表中 p 为增加后的处理器数量，n 为扩大后的数据集大小；各项均按每处理器计。"} 计算量随 \(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 个处理器的执行时间为： $$ \begin{aligned} T_{32}&=\text{计算时间}+\text{通信时间}\\ &=\frac{0.8T_4}{8}+(0.2T_4)\sqrt8\\ &=0.1T_4+0.57T_4=0.67T_4 \end{aligned} $$ 于是加速比为： $$ \text{加速比}=\frac{T_4}{T_{32}}=\frac{T_4}{0.67T_4}=1.49 $$ 通信时间占比从 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 倍的结论不受此处方向错误影响。