并行系统中稀疏矩阵计算的非结构化网格排序策略

### 稀疏矩阵计算中无结构网格排序策略研究 #### 1. 引言 计算机借助数学模型解决复杂问题和模拟复杂过程的能力，使其成为现代科学与工程领域不可或缺的工具。在大规模现实应用的计算机模拟中，通常需要在有限区域内求解一组非线性偏微分方程（PDEs）。例如，美国能源部重大挑战项目的一个重点领域是设计未来的加速器，如散裂中子源（SNS）。斯坦福直线加速器中心（SLAC）的研究人员需要对具有大长宽比的复杂射频四极（RFQ）腔进行建模。 目前，无结构网格被用于解决大型计算域中的小特征问题，未来还将添加动态网格自适应技术以提高效率。电磁学中的PDEs通过有限元法（FEM）进行离散化，从而得到广义特征值问题 $Kx = \lambda Mx$，其中 $K$ 和 $M$ 分别是刚度矩阵和质量矩阵，且均为稀疏矩阵。在典型的腔模型中，自由度数量约为一百万。对于如此大规模的特征值问题，直接求解技术很快会达到内存限制，因此最广泛使用的方法是Krylov子空间方法，如Lanczos或Jacobi - Davidson方法。在所有基于Krylov的算法中，稀疏矩阵向量乘法（SPMV）必须反复执行，因此SPMV的效率通常决定了特征求解器的速度。SPMV也是大规模数值模拟中使用最频繁的核心操作之一。 在单处理器机器上，求解此类复杂现实问题的数值解可能极其耗时，这推动了越来越强大的并行（多处理器）超级计算机的发展。然而，许多值得模拟的系统具有无结构和动态的特性，这使得它们的高效并行实现成为一项艰巨的任务。此外，基于深层内存层次结构的现代计算机体系结构，只有在用户关注数据的正确分布和放置时才能表现出可接受的性能。单处理器性能关键取决于局部性的利用，而如果数据分区不当导致过多的通信和/或数据迁移，并行性能将显著下降。传统方法是使用像MeTiS这样的分区工具，并对得到的分区进行枚举策略后处理以增强局部性。尽管分区和局部性优化可以被视为两个独立的问题，但实际应用中两者之间往往存在复杂的相互作用。 #### 2. 分区与线性化 空间填充曲线已被证明是解决N体和FEM模拟中某些问题的一种优雅而统一的线性化方法。将高维空间结构线性化，即将其映射到一维超空间，有两种应用方式：一是其“保持局部性”的特性与给定的内存层次结构完美匹配；二是对连续线性对象进行分区非常简单。在实验中，我们对这两种策略进行了一些修改。以下简要介绍所使用的两种枚举技术和通用图分区器。 - **Cuthill - McKee算法（CM）**：在FEM离散化中，顶点的特定枚举在很大程度上控制着所得刚度矩阵的稀疏模式。矩阵的带宽或轮廓对线性系统和特征求解器的效率有显著影响。Cuthill和McKee提出了一种基于图论思想的简单算法。从度最小的顶点开始，首先构建与该顶点“距离”逐渐增加的层次，然后按层次进行枚举，每个层次内按顶点度递增的顺序进行。该方法有多种变体，最流行的是反向Cuthill - McKee（RCM），它从最后一层中度最小的顶点重新开始层次构建。CM算法类实现相对简单，并且主要在纯图结构上运行，即底层图不一定来自三角形网格。 - **自回避行走（SAW）**：这是最近提出的一种基于网格（而非基于几何）的技术，与空间填充曲线具有相似的应用领域。三角形网格上的SAW是对三角形的一种枚举，使得SAW中连续的两个三角形共享一条边或一个顶点，即SAW中没有跳跃。可以证明，在任意无结构网格上存在具有更特殊性质的行走，并且有一个复杂度与网格中三角形数量呈线性关系的算法来构造它们。此外，SAW适用于分层粗化和细化，即只需在网格自适应发生的区域重建，因此可以很容易地并行化。与CM不同，SAW不是专门为顶点枚举设计的技术，因此它不能在三角形网格的裸图结构上运行，这意味着SAW的构造成本较高，但可以从给定的SAW导出几种不同的顶点枚举。 - **图分区（MeTiS）**：在过去5 - 10年中，已经开发并实现了一些优秀的并行图分区算法，这些算法速度极快，同时能提供良好的负载平衡质量和低边割。其中最流行的可能是MeTiS，它属于多级分区器类。MeTiS使用重边匹配方案合并顶点和边来减小图的大小，对最粗化的图应用贪心图生长算法进行分区，然后使用边界贪心和Kernighan - Lin细化的组合进行反粗化，以构建原始图的分区。 #### 3. 实验结果 为了执行SPMV（$y \leftarrow Ax$），假设矩阵 $A$ 的非零元素以压缩稀疏行格式存储，密集向量 $x$ 以单位步长顺序存储在内存中。网格元素/顶点的不同编号会导致 $A$ 的非零模式不同，进而导致访问 $x$ 元素的模式不同。此外，在分布式内存机器上，这意味着不同的通信量。 实验测试网格是由Triangle生成的二维Delaunay三角剖分，形状像字母“A”，包含661,054个顶点和1,313,099个三角形。底层矩阵是通过为网格中边的顶点端点 $(v_1, v_2)$ 对应的每个（行，列）条目分配一个随机值来组装的，这模拟了每个顶点需要与其最近邻通信的模板计算。最终矩阵极其稀疏，仅包含2,635,207个非零元素。SPMV所需的浮点运算次数是非零元素数量的两倍（对于测试矩阵为5,270,414次）。 ##### 3.1 分布式内存实现 在实验中，使用基于MPI实现的Aztec中的并行SPMV例程。矩阵 $A$ 被划分为行块，每个块分配给一个处理器。有两个特别感兴趣的例程：$AZ\_transform$ 和 $AZ\_matvec\_mult$。前者初始化数据结构和通信调度，后者执行矩阵向量乘法。表1报告了这些例程在NERSC的450 MHz Cray T3E上的运行时间。原始自然排序（ORIG）最慢，在分布式内存机器上显然不可接受。对于关键内核例程 $AZ\_matvec\_mult$，RCM略快于S