MX Direct Solver——高性能稀疏线性求解器

2023-03-21 14:51 作者:Mack丨迈曦科技 0人读过 | 我要投稿

MXDSV是基于直接方法的大规模稀疏线性方程求解器，目前支持非对称、结构对称以及对称方程组，实数与复数方程组等求解功能。

求解形式为Ax=b的线性方程组的直接方法是基于高斯消元计算A=LU，其中A和U分别为下三角和上三角矩阵，计算系数矩阵A的三角分解也称为LU分解。分解后，通过求解三角形系统 Ly=b 和Ux=y即可简单求解原系统。如果A是对称的，则通过Cholesky分解来计算形式A=LLT或A=LDLT的分解，其中A是一个下三角的矩阵(在A=LDLT分解的情况下，为单位下三角矩阵)，并且D是一个对角矩阵。

图1和图2分别展示了稠密矩阵的LU分解和Cholesky分解的一组常用公式。在此基础上，通过在这些算法中重新排列循环，可以得到其他数学上等价的公式。这些算法同样适用于稀疏矩阵—其大部分元素是零。

例如，如果A[j,i]除法步骤中的分子为零，则不需要执行该操作。类似地，如果A[j,i]或A[i,k](如果是对称的，则A[k,i])为零，则可以省略更新步骤。

图1 列优先的密集矩阵LU分解算法

图2 列优先的密集矩阵Cholesky分解算法

当A是稀疏的时，三角矩阵L和U通常会在额外的位置产生非零元，这种现象被称为填充（fill-in），这会导致使用直接方法求解稀疏系统的内存和时间需求相对于系统产生超线性增长。尽管直接方法需要很高的内存，但由于其通用性和鲁棒性，在许多工程实际应用中仍然是首选方法。在需要求解多个右端向量并且系数矩阵不变的应用中，直接方法更是求解者的选择，因为一次数值分解的开销可以平摊在多个快速的三角矩阵求解上。

01 算法线路

稀疏线性系统的直接求解过程通常包括四个阶段。两个计算阶段，因式分解和三角求解已经提到过。有些时候，因子中非零元的数目是系数矩阵的行和列的初始排列的函数。在许多并行的稀疏分解形式中，这种排列也会对负载平衡产生影响。因此，直接求解稀疏线性系统的第一步是应用启发式算法来计算矩阵的理想排列，这一步被称为排序。稀疏矩阵可以看作是图的邻接矩阵，排序启发式通常使用矩阵的图形视图，并按照特定的顺序标记顶点，这相当于计算具有理想属性的系数矩阵的排列。第二个阶段，称为符号分解，计算因子的非零元模式（分布）。

在实际计算因子之前了解因子的非零模式是有用的，原因有几个：

(a)在符号分解过程中，可以预测数值分解的内存需求；

(b)由于预先知道非零的数量和位置，在数值分解过程中可以避免大量的间接寻址，从而提高性能；

(c)在并行实现中，符号分解有助于数据和计算在处理单元之间的分配。

排序和符号分解阶段也称为预处理或分析步骤。迭代细化的第五个阶段有时在求解阶段之后使用，以提高求解的精度。

在实际计算中，数值分解通常消耗最多的内存和时间。许多应用涉及到分解几个具有不同数值但具有相同稀疏结构的矩阵，在这种情况下，排序和符号分解步骤的部分或全部结果都可以重用。这对于并行稀疏求解器也是有利的，因为并行排序和符号分解通常是可扩展性较差的。将这些步骤的成本分摊到几个分解步骤中，有助于保持求解器的整体可伸缩性接近于数值分解的可伸缩性。三角求解的并行化高度依赖于数值分解阶段的并行化。数值因式分解的并行形式规定了因子如何在并行任务中分配。随后的三角求解步骤必须使用一个并行化方案来处理这种数据分布，特别是在分布式内存并行环境中。

02 异构计算

目前对于线性方程直接法中数值分解的并行加速工作主要基于两个方面：一个是节点并行，另一个是树并行。很多直接求解器使用GPU加速数值分解都是基于节点并行的方式，也就是说利用GPU对部分dense kernel（blas or lapack）进行加速计算。在许多问题中这种辅助计算的方法并不能充分利用GPU的计算资源，所以计算效率难以得到大幅提升。然而，GPU拥有大量可并发的线程和高访存带宽，这为我们提供一种新的并行策略：同时使用节点并行和树并行。因此，我们提出了一种面向多GPU的多级并行策略对稀疏分解进行加速计算，并在此基础上针对有限元分析设计高效地求解算法。

03MXDSV性能表现

我们从SuiteSparse Matrix Collection (http://sparse.tamu.edu/)矩阵集中选取了前10个较大的稀疏矩阵作为测试集，与知名商业求解器MKL PARDISO进行求解性能对比。