大规模稀疏矩阵PARDISO求解方法介绍.pdf

高性能计算 13 大规模稀疏矩阵PARDISO求解方法介绍 于 超 英特尔亚太研发中心 200241上海 chao.yu@ 摘要： 大规模稀疏矩阵的求解是高性能计算中的一个常见问题。本文介绍了用直接法（Direct ® Sparse Solver）求解矩阵的一些问题以及使用IntelMKL PARDISO接口求解稀疏矩阵的方法。 1.引言 但是对于稀疏矩阵来说，在求解的过程中会遇 大规模稀疏方程的求解是工程计算中常见的问 到两个问题：非零元素填入 （Fill-in）与矩阵行列的 题, 一些典型应用包括有限元求解,积分方程，求矩 重排 （Reordering）。我们用一个具体的例子来说明 阵特征值，最优化问题等。这些问题要求稀疏方程 这个问题: 的求解方法具有如下一些特性： （1）求解方法性能 假设有如下的稀疏矩阵方程：Ax = b。A为对称 较高。稀疏矩阵维数相对较大,对求解性能要求高； 正定稀疏矩阵: （2）高效的矩阵存储方法存储稀疏矩阵与中间计算 结果； （3）求解过程具有稳定性,特别对一些病态的 矩阵能够得到正确结果。 [3] PARDISO是在共享内存机器上实现的稀疏 矩阵的求解方法，对于一些大规模的计算问题, PARDISO的算法表现了非常好的计算效率与并行 性。一些数值测试表明，随着计算节点数目增加, A(*) 说明A中相应的元素0. 矩阵A是一个绝大多 PARDISO具有接近线性的加速比例[3]。下面我们结 数元素都为零的稀疏矩阵. 其相应的LU分解(对称正定 合Intel MKL函数库PARDISO 的接口，介绍直接法来 T 矩阵,称为Cholesky分解)为 A =LL ： 求解矩阵的PARDISO方法。 2. DirectSparse Solver 的求解方法 用直接方法来求解矩阵，一个最基本的方法是 将矩阵分解为上下三角矩阵。也就是说，对于线性 方程： Ax =b 我们需要找到一个下三角的矩阵L与上三角的矩 从分解的结果中，我们可以看出，尽管矩阵 A 阵U(通常称为LU分解), 使得 为稀疏矩阵，但是其分解后上下三角矩阵L却有较多 A =LU -LUx=b 的非零元素.。这样，如果直接计算L来求解方程，我 这样，我们只需两步求解上下三角矩阵来求解 们实际会与稠密矩阵相似的计算复杂度与数据存储 原有方程： 量。 1.求解 Ly = b. 如果矩阵A中的某一零元素，对其分解后，矩阵 2.求解Ux = y. L相应位置产生非零元素。我们称之为非零元素填入 求解 Ly = b或Ux = y可以用简单的正向与反向的 (Fill-in)。从计算的角度，一种有效的方法是我们去 方法求解上下三角矩阵。 遍历A非零元素，同时尽量减少分解后矩阵L中非零 14 《高性能计算发展与应用》 2006年第四期 总第十七期 元素的个数，这样，我们求解时可以只需要考虑L中