高性能计算
13
大规模稀疏矩阵ＰＡＲＤＩＳＯ求解方法介绍
于 超
英特尔亚太研发中心 上海 ２００２４１ ｃｈａｏ．ｙｕ＠ｉｎｔｅｌ．ｃｏｍ
摘要：
大规模稀疏矩阵的求解是高性能计算中的一个常见问题。本文介绍了用直接法（Ｄｉｒｅｃｔ Ｓｐａｒｓｅ Ｓｏｌｖｅｒ）求解矩阵的一些问题以及使用ＩｎｔｅｌＭＫＬ ＰＡＲＤＩＳＯ接口求解稀疏矩阵的方法。
１． 引言
大规模稀疏方程的求解是工程计算中常见的问题， 一些典型应用包括有限元求解， 积分方程，求矩阵特征值，最优化问题等。这些问题要求稀疏方程的求解方法具有如下一些特性：（１）求解方法性能较高。稀疏矩阵维数相对较大，对求解性能要求高；（２）高效的矩阵存储方法存储稀疏矩阵与中间计算结果；（３）求解过程具有稳定性，特别对一些病态的矩阵能够得到正确结果。
ＰＡＲＤＩＳＯ［３］是在共享内存机器上实现的稀疏矩阵的求解方法，对于一些大规模的计算问题， ＰＡＲＤＩＳＯ的算法表现了非常好的计算效率与并行性。一些数值测试表明，随着计算节点数目增加， ＰＡＲＤＩＳＯ具有接近线性的加速比例［３］。下面我们结合Ｉｎｔｅｌ ＭＫＬ函数库ＰＡＲＤＩＳＯ 的接口，介绍直接法来求解矩阵的ＰＡＲＤＩＳＯ方法。
但是对于稀疏矩阵来说，在求解的过程中会遇到两个问题：非零元素填入（Ｆｉｌｌ－ｉｎ）与矩阵行列的重排（Ｒｅｏｒｄｅｒｉｎｇ）。我们用一个具体的例子来说明这个问题：
假设有如下的稀疏矩阵方程：Ａｘ ＝ ｂ。Ａ 为对称正定稀疏矩阵：
Ａ（＊） 说明Ａ中相应的元素０． 矩阵Ａ是一个绝大多数元素都为零的稀疏矩阵． 其相应的ＬＵ分解（对称正定矩阵，称为Ｃｈｏｌｅｓｋｙ分解）为 Ａ ＝ LL T ：
２． Ｄｉｒｅｃｔ Ｓｐａｒｓｅ Ｓｏｌｖｅｒ 的求解方法
用直接方法来求解矩阵，一个最基本的方法是将矩阵分解为上下三角矩阵。也就是说，对于线性方程：
从分解的结果中，我们可以看出，尽管矩阵 Ａ 为稀疏矩阵，但是其分解后上下三角矩阵Ｌ却有较多的非零元素．。这样，如果直接计算Ｌ来求解方程，我们实际会与稠密矩阵相似的计算复杂度与数据存储量。
如果矩阵Ａ中的某一零元素，对其分解后，矩阵 Ｌ相应位置产生非零元素。我们称之为非零元素填入（Ｆｉｌｌ－ｉｎ）。从计算的角度，一种有效的方法是我们去遍历Ａ非零元素，同时尽量减少分解后矩阵Ｌ中非零
Ａｘ ＝ ｂ
我们需要找到一个下三角的矩阵Ｌ与上三角的矩阵Ｕ（通常称为ＬＵ分解）， 使得
Ａ ＝ ＬＵ －＞ ＬＵｘ＝ｂ
这样，我们只需两步求解上下三角矩阵来求解原有方程：
１． 求解 Ｌｙ ＝ ｂ．２． 求解 Ｕｘ ＝ ｙ．
求解 Ｌｙ ＝ ｂ 或 Ｕｘ ＝ ｙ 可以用简单的正向与反向的方法求解上下三角矩阵。
少量非零元素。
在实际计算中，一种常用的方式就是对Ａ的行与列进行交换。我们可以用一个交换矩阵 Ｐ 来表示矩阵的行列的交换。例如，在上面的例子中，我们如果交换矩阵Ａ的第一行与第五行以及第一列与第五列，得到矩阵 Ｂ ＝
：
解方程，这样方程求解包括如下求解步骤：
１． 求交换矩阵Ｐ： Ｂ ＝ ＰＡＰＴ２． 分解矩阵 Ｂ ＝ ＬＵ３． 求解 Ｌｙ ＝ Ｐｂ４． 求解 Ｕｚ ＝ ｙ５． 计算 Ｘ ＝ ＰＴｚ
３． ＰＡＲＤＩＳＯ 求解方法
ＰＡＲＤＩＳＯ求解过程
ＰＡＲＤＩＳＯ为Ｂａｓｅｌ大学提供的一个稀疏矩阵接口。 Ｉｎｔｅｌ Ｍａｔｈ Ｋｅｒｎｅｌ Ｌｉｂｒａｒｙ提供优化版本。根据如
Ｂ仅仅是Ａ交换行列后得到的矩阵，因此Ｂ与Ａ有相同的非零元素，但是，我们对矩阵Ｂ进行分解得到 Ｂ ＝ LL T ：
上讨论的求解稀疏矩阵的方法，ＰＡＲＤＩＳＯ对应求解过程包括如下步骤：
１． 矩阵重排与符号分解（Ｒｅｏｒｄｅｒｉｎｇ ａｎｄ Ｓｙｍｂｏｌｉｃ Ｆａｃｔｏｒｉｚａｔｉｏｎ）：ＰＡＲＤＩＳＯ Ｓｏｌｖｅｒ根据不同的矩阵类型，计算不同类型的行列交换矩阵Ｐ与对角矩阵Ｄ，对Ａ矩阵进行交换重排。新得到的矩阵括尽量少的非零元素。
２． 矩阵ＬＵ 分解： 对
可以看出。Ｂ分解后的Ｌ矩阵的非零元素比Ａ分解后的矩阵要少很多，也就说，Ｂ分解的非零元素填入要比Ａ少很多。进而矩阵Ｂ的存储与计算的复杂度要比Ａ分解要明显减少。所以一个高效稀疏矩阵的求解首先须找出一个或多个交换矩阵
，Ｐ能够减少非零元素填
入。
如上的方法是基于一个对称的正定矩阵。对于一般化稀疏矩阵，我们有类似的工作。首先用交换矩阵Ｐ对原有Ａ进行交换，然后进行ＬＵ分解，后再求
进行 ＬＵ 分解。
３． 方程求解与迭代：根据ＬＵ分解的结果，求解方程，如果对结果的精度有进一步要求，使用迭代法进一步提高解精度。
４． 迭代结束，释放计算过程的内存。
除ＰＡＲＤＩＳＯ外，Ｉｎｔｅｌ ＭＫＬ 还提供一个较为容易使用的Ｄｉｒｅｃｔ Ｓｐａｒｓｅ Ｓｏｌｖｅｒ（ＤＳＳ）的接口。
矩阵类型
ＰＡＲＤＩＳＯ函数的接口对多种的矩阵类型提供支持，包括了实、复数、对称或不对称的矩阵。具体可以如图１所示：
分解后会包
图１ ＰＡＲＤＩＳＯ支持的矩阵类型
高性能计算
15
矩阵存储
目前，稀疏矩阵存在多种存储格式。ＰＡＲＤＩＳＯ接口中使用以行为主的存储方法。该方法以行为单位存储每个非零数据。对于一个稀疏矩阵Ａ，ＰＡＲＤＩＳＯ对矩阵的存储包括了三个数组：
• • •
ｖａｌｕｅｓ – 矩阵Ａ的实数或复数非零数据． Ａ的非ｃｏｌｕｍｎｓ － ｖａｌｕｅｓ 中每个元素所在矩阵的列． ｒｏｗｉｎｄｅｘ– 给出每一行的元素在ｖａｌｕｅ 中的位置．
零数据通过下面ｃｏｌｕｍｎｓ与ｒｏｗｉｎｄｅｘ映射到ｖａｌｕｅｓ 数组中．
性能数据
一些数值测试表明，ＰＡＲＤＩＳＯ是目前最快的线性稀疏矩阵的求解方法之一［２］。下图２给出ＰＡＲＤＩＳＯ与其他稀疏矩阵求解法的性能比较。图中数据是对一些较大测试矩阵，其他求解法时间与ＰＡＲＤＩＳＯ的时间比（其他求解器的求解时间／ ＰＡＲＳＤＩＳＯ求解的时间）。有关这方面的更多数据，
可参见文［２］的测试。
图２ ＰＡＲＤＩＳＯ与其他求解器的性能对比
４． 小结
本文讨论了一个广为使用的稀疏矩阵的求解接口ＰＡＲＤＩＳＯ。Ｉｎｔｅｌ ＭＫＬ 函数提供了高效的ＰＡＲＤＩＳＯ求解方法实现，此外，Ｉｎｔｅｌ ＭＫＬ 函数库中还包括其他高性能的优化函数，如向量与矩阵运算的ＢＬＡＳ、
ＬＡＰＡＣＫ函数。多维傅立叶变换函数，向量数学函数（ＶＭＬ）以及随机数生产函数等等。进一步的内容可从如下网站获得：ｈｔｔｐ：／／ｗｗｗ．ｉｎｔｅｌ．ｃｏｍ／ｓｏｆｔｗａｒｅ／ｐｒｏｄｕｃｔｓ／ＭＫＬ
参考文献
［１］ Ｉｎｔｅｌ Ｍａｔｈ Ｋｅｒｎｅｌ Ｌｉｂｒａｒｙ： ｈｔｔｐ：／／ｗｗｗ．ｉｎｔｅｌ．ｃｏｍ／ｓｏｆｔｗａｒｅ／ｐｒｏｄｕｃｔｓ／ＭＫＬ
［２］ Ｎ Ｉ Ｍ Ｇｏｕｌｄ， Ｙ Ｈｕ． Ｊ Ａ Ｓｃｏｔｔ． Ｃｏｍｐｌｅｔｅ ｒｅｓｕｌｔｓ ｆｒｏｍ ａ ｎｕｍｅｒｉｃａｌ ｅｖａｌｕａｔｉｏｎ ｏｆ ｓｐａｒｓｅ ｄｉｒｅｃｔ ｓｏｌｖｅｒｓ ｆｏｒ ｔｈｅ ｓｏｌｕｔｉｏｎ ｏｆ ｌａｒｇｅ， ｓｐａｒｓｅ， ｓｙｍｍｅｔｒｉｃ ｌｉｎｅａｒ ｓｙｓｔｅｍｓ ｏｆ ｅｑｕａｔｉｏｎｓ． ｈｔｔｐ：／／ｗｗｗ．ｎｕｍｅｒｉｃａｌ．ｒｌ．ａｃ．ｕｋ／ｒｅｐｏｒｔｓ／ｒｅｐｏｒｔｓ．ｓｈｔｍｌ［３］ Ｓｏｆｔｗａｒｅ ＰＡＲＤＩＳＯ： ｈｔｔｐ：／／ｗｗｗ．ｃｏｍｐｕｔａｔｉｏｎａｌ．ｕｎｉｂａｓ．ｃｈ／ｃｓ／ｓｃｉｃｏｍｐ／ｓｏｆｔｗａｒｅ／ｐａｒｄｉｓｏ／
高性能计算
13
大规模稀疏矩阵ＰＡＲＤＩＳＯ求解方法介绍
于 超
英特尔亚太研发中心 上海 ２００２４１ ｃｈａｏ．ｙｕ＠ｉｎｔｅｌ．ｃｏｍ
摘要：
大规模稀疏矩阵的求解是高性能计算中的一个常见问题。本文介绍了用直接法（Ｄｉｒｅｃｔ Ｓｐａｒｓｅ Ｓｏｌｖｅｒ）求解矩阵的一些问题以及使用ＩｎｔｅｌＭＫＬ ＰＡＲＤＩＳＯ接口求解稀疏矩阵的方法。
１． 引言
大规模稀疏方程的求解是工程计算中常见的问题， 一些典型应用包括有限元求解， 积分方程，求矩阵特征值，最优化问题等。这些问题要求稀疏方程的求解方法具有如下一些特性：（１）求解方法性能较高。稀疏矩阵维数相对较大，对求解性能要求高；（２）高效的矩阵存储方法存储稀疏矩阵与中间计算结果；（３）求解过程具有稳定性，特别对一些病态的矩阵能够得到正确结果。
ＰＡＲＤＩＳＯ［３］是在共享内存机器上实现的稀疏矩阵的求解方法，对于一些大规模的计算问题， ＰＡＲＤＩＳＯ的算法表现了非常好的计算效率与并行性。一些数值测试表明，随着计算节点数目增加， ＰＡＲＤＩＳＯ具有接近线性的加速比例［３］。下面我们结合Ｉｎｔｅｌ ＭＫＬ函数库ＰＡＲＤＩＳＯ 的接口，介绍直接法来求解矩阵的ＰＡＲＤＩＳＯ方法。
但是对于稀疏矩阵来说，在求解的过程中会遇到两个问题：非零元素填入（Ｆｉｌｌ－ｉｎ）与矩阵行列的重排（Ｒｅｏｒｄｅｒｉｎｇ）。我们用一个具体的例子来说明这个问题：
假设有如下的稀疏矩阵方程：Ａｘ ＝ ｂ。Ａ 为对称正定稀疏矩阵：
Ａ（＊） 说明Ａ中相应的元素０． 矩阵Ａ是一个绝大多数元素都为零的稀疏矩阵． 其相应的ＬＵ分解（对称正定矩阵，称为Ｃｈｏｌｅｓｋｙ分解）为 Ａ ＝ LL T ：
２． Ｄｉｒｅｃｔ Ｓｐａｒｓｅ Ｓｏｌｖｅｒ 的求解方法
用直接方法来求解矩阵，一个最基本的方法是将矩阵分解为上下三角矩阵。也就是说，对于线性方程：
从分解的结果中，我们可以看出，尽管矩阵 Ａ 为稀疏矩阵，但是其分解后上下三角矩阵Ｌ却有较多的非零元素．。这样，如果直接计算Ｌ来求解方程，我们实际会与稠密矩阵相似的计算复杂度与数据存储量。
如果矩阵Ａ中的某一零元素，对其分解后，矩阵 Ｌ相应位置产生非零元素。我们称之为非零元素填入（Ｆｉｌｌ－ｉｎ）。从计算的角度，一种有效的方法是我们去遍历Ａ非零元素，同时尽量减少分解后矩阵Ｌ中非零
Ａｘ ＝ ｂ
我们需要找到一个下三角的矩阵Ｌ与上三角的矩阵Ｕ（通常称为ＬＵ分解）， 使得
Ａ ＝ ＬＵ －＞ ＬＵｘ＝ｂ
这样，我们只需两步求解上下三角矩阵来求解原有方程：
１． 求解 Ｌｙ ＝ ｂ．２． 求解 Ｕｘ ＝ ｙ．
求解 Ｌｙ ＝ ｂ 或 Ｕｘ ＝ ｙ 可以用简单的正向与反向的方法求解上下三角矩阵。
少量非零元素。
在实际计算中，一种常用的方式就是对Ａ的行与列进行交换。我们可以用一个交换矩阵 Ｐ 来表示矩阵的行列的交换。例如，在上面的例子中，我们如果交换矩阵Ａ的第一行与第五行以及第一列与第五列，得到矩阵 Ｂ ＝
：
解方程，这样方程求解包括如下求解步骤：
１． 求交换矩阵Ｐ： Ｂ ＝ ＰＡＰＴ２． 分解矩阵 Ｂ ＝ ＬＵ３． 求解 Ｌｙ ＝ Ｐｂ４． 求解 Ｕｚ ＝ ｙ５． 计算 Ｘ ＝ ＰＴｚ
３． ＰＡＲＤＩＳＯ 求解方法
ＰＡＲＤＩＳＯ求解过程
ＰＡＲＤＩＳＯ为Ｂａｓｅｌ大学提供的一个稀疏矩阵接口。 Ｉｎｔｅｌ Ｍａｔｈ Ｋｅｒｎｅｌ Ｌｉｂｒａｒｙ提供优化版本。根据如
Ｂ仅仅是Ａ交换行列后得到的矩阵，因此Ｂ与Ａ有相同的非零元素，但是，我们对矩阵Ｂ进行分解得到 Ｂ ＝ LL T ：
上讨论的求解稀疏矩阵的方法，ＰＡＲＤＩＳＯ对应求解过程包括如下步骤：
１． 矩阵重排与符号分解（Ｒｅｏｒｄｅｒｉｎｇ ａｎｄ Ｓｙｍｂｏｌｉｃ Ｆａｃｔｏｒｉｚａｔｉｏｎ）：ＰＡＲＤＩＳＯ Ｓｏｌｖｅｒ根据不同的矩阵类型，计算不同类型的行列交换矩阵Ｐ与对角矩阵Ｄ，对Ａ矩阵进行交换重排。新得到的矩阵括尽量少的非零元素。
２． 矩阵ＬＵ 分解： 对
可以看出。Ｂ分解后的Ｌ矩阵的非零元素比Ａ分解后的矩阵要少很多，也就说，Ｂ分解的非零元素填入要比Ａ少很多。进而矩阵Ｂ的存储与计算的复杂度要比Ａ分解要明显减少。所以一个高效稀疏矩阵的求解首先须找出一个或多个交换矩阵
，Ｐ能够减少非零元素填
入。
如上的方法是基于一个对称的正定矩阵。对于一般化稀疏矩阵，我们有类似的工作。首先用交换矩阵Ｐ对原有Ａ进行交换，然后进行ＬＵ分解，后再求
进行 ＬＵ 分解。
３． 方程求解与迭代：根据ＬＵ分解的结果，求解方程，如果对结果的精度有进一步要求，使用迭代法进一步提高解精度。
４． 迭代结束，释放计算过程的内存。
除ＰＡＲＤＩＳＯ外，Ｉｎｔｅｌ ＭＫＬ 还提供一个较为容易使用的Ｄｉｒｅｃｔ Ｓｐａｒｓｅ Ｓｏｌｖｅｒ（ＤＳＳ）的接口。
矩阵类型
ＰＡＲＤＩＳＯ函数的接口对多种的矩阵类型提供支持，包括了实、复数、对称或不对称的矩阵。具体可以如图１所示：
分解后会包
图１ ＰＡＲＤＩＳＯ支持的矩阵类型
高性能计算
15
矩阵存储
目前，稀疏矩阵存在多种存储格式。ＰＡＲＤＩＳＯ接口中使用以行为主的存储方法。该方法以行为单位存储每个非零数据。对于一个稀疏矩阵Ａ，ＰＡＲＤＩＳＯ对矩阵的存储包括了三个数组：
• • •
ｖａｌｕｅｓ – 矩阵Ａ的实数或复数非零数据． Ａ的非ｃｏｌｕｍｎｓ － ｖａｌｕｅｓ 中每个元素所在矩阵的列． ｒｏｗｉｎｄｅｘ– 给出每一行的元素在ｖａｌｕｅ 中的位置．
零数据通过下面ｃｏｌｕｍｎｓ与ｒｏｗｉｎｄｅｘ映射到ｖａｌｕｅｓ 数组中．
性能数据
一些数值测试表明，ＰＡＲＤＩＳＯ是目前最快的线性稀疏矩阵的求解方法之一［２］。下图２给出ＰＡＲＤＩＳＯ与其他稀疏矩阵求解法的性能比较。图中数据是对一些较大测试矩阵，其他求解法时间与ＰＡＲＤＩＳＯ的时间比（其他求解器的求解时间／ ＰＡＲＳＤＩＳＯ求解的时间）。有关这方面的更多数据，
可参见文［２］的测试。
图２ ＰＡＲＤＩＳＯ与其他求解器的性能对比
４． 小结
本文讨论了一个广为使用的稀疏矩阵的求解接口ＰＡＲＤＩＳＯ。Ｉｎｔｅｌ ＭＫＬ 函数提供了高效的ＰＡＲＤＩＳＯ求解方法实现，此外，Ｉｎｔｅｌ ＭＫＬ 函数库中还包括其他高性能的优化函数，如向量与矩阵运算的ＢＬＡＳ、
ＬＡＰＡＣＫ函数。多维傅立叶变换函数，向量数学函数（ＶＭＬ）以及随机数生产函数等等。进一步的内容可从如下网站获得：ｈｔｔｐ：／／ｗｗｗ．ｉｎｔｅｌ．ｃｏｍ／ｓｏｆｔｗａｒｅ／ｐｒｏｄｕｃｔｓ／ＭＫＬ
参考文献
［１］ Ｉｎｔｅｌ Ｍａｔｈ Ｋｅｒｎｅｌ Ｌｉｂｒａｒｙ： ｈｔｔｐ：／／ｗｗｗ．ｉｎｔｅｌ．ｃｏｍ／ｓｏｆｔｗａｒｅ／ｐｒｏｄｕｃｔｓ／ＭＫＬ
［２］ Ｎ Ｉ Ｍ Ｇｏｕｌｄ， Ｙ Ｈｕ． Ｊ Ａ Ｓｃｏｔｔ． Ｃｏｍｐｌｅｔｅ ｒｅｓｕｌｔｓ ｆｒｏｍ ａ ｎｕｍｅｒｉｃａｌ ｅｖａｌｕａｔｉｏｎ ｏｆ ｓｐａｒｓｅ ｄｉｒｅｃｔ ｓｏｌｖｅｒｓ ｆｏｒ ｔｈｅ ｓｏｌｕｔｉｏｎ ｏｆ ｌａｒｇｅ， ｓｐａｒｓｅ， ｓｙｍｍｅｔｒｉｃ ｌｉｎｅａｒ ｓｙｓｔｅｍｓ ｏｆ ｅｑｕａｔｉｏｎｓ． ｈｔｔｐ：／／ｗｗｗ．ｎｕｍｅｒｉｃａｌ．ｒｌ．ａｃ．ｕｋ／ｒｅｐｏｒｔｓ／ｒｅｐｏｒｔｓ．ｓｈｔｍｌ［３］ Ｓｏｆｔｗａｒｅ ＰＡＲＤＩＳＯ： ｈｔｔｐ：／／ｗｗｗ．ｃｏｍｐｕｔａｔｉｏｎａｌ．ｕｎｉｂａｓ．ｃｈ／ｃｓ／ｓｃｉｃｏｍｐ／ｓｏｆｔｗａｒｅ／ｐａｒｄｉｓｏ／