转载：Intel MKL 稀疏矩阵求解PARDISO 函数

Posted 2020-10-01 Kejia Pan

Intel MKL提供了针对稀疏矩阵求解的PARDISO 接口，它是在共享内存机器上，实现的稀疏矩阵的直接求解方法，对于一些大规模的计算问题, PARDISO的算法表现了非常好的计算效率与并行性。一些数值测试表明，随着计算节点数目增加, PARDISO具有接近线性的加速比例。

 PARDISO对应求解过程包括如下步骤：

1. 矩阵重排与符号分解(Reordering and Symbolic Factorization)：PARDISO Solver根据不同的矩阵类型，计算不同类型的行列交换矩阵P与对角矩阵D，对A矩阵进行交换重排。新得到的矩阵分解后会包括尽量少的非零元素。

2. 矩阵LU 分解: 对 进行 LU 分解。

3. 方程求解与迭代：根据LU分解的结果，求解方程，如果对结果的精度有进一步要求，使用迭代法进一步提高解精度。

4. 迭代结束，释放计算过程的内存。

使用PARDISO 的时候，可能会有一些常见的问题：

第一，Paridso 提示内存不足：

  出现这类问题的时候，可以首先检查一下Pardiso 对求解该问题的内存需求，Paridiso计算时，可以通过下面的数据求得：

   max(iparm(15), iparm(16)+iparm(17))

 可以对比一下这个数据，查看系统的内存是否满足需求。

Paridso 同时支持，in-core 与 out-of-core的计算。 如果， 计算的数据太大，而不能完全在内存求解的时候，可以的使用out-of-core 的pardiso（设置 iparm(60) 参数）。 Out-of-core 的计算会将中间计算数据保存于硬盘上，从而能够解决一些大的计算问题。

实际中，还常常遇到的一个问题是，许多应用是32位的程序， 这样，即使使用out-of-core 的pardiso来求解，仍然会受到32位的地址空间的限制。 如果计算数据非常的大，需要改写为64 位的计算程序。

第二，检查输入数据的合法性：

使用Pardiso 在进行计算的时候，常常会出现中间计算错误。由于Pardiso采用CSR格式的压缩存储的矩阵。 很多情况下，计算错误是由于输入了不合法的计算数据而导致。对于这类问题，可以在调用Pardiso 的时候， 进行输入数据的检查（设置iparm(27) 的参数）， Paridso如果发现输入数据的错误，会给出错误提示。这类检查，可以帮助发现一些简单的,特别是与输入数据的索引相关的输入错误。

第三，使用缺省参数：

Pardiso 中提供了丰富的输入参数选项。 用户在调用的时候，需要确保正确的输入参数。 很多在计算过程中发生的错误，往往与不正确的输入参数相关。 一个常用的检查方法是输入缺省的paridso 的参数（iparm(1) =0），Paridso 使用缺省参数进行计算，来检验程序的正确性。

第四，在C/C++语言的调用Pardiso：

在Intel MKL 函数手册中， Pardiso 相关参数的说明是以Fortran语言的形式给出。 如果我们在C/C++语言中，调用Pardiso 函数，需要注意输入数据的数组下标。 C语言中对应的数组下标是从0开始，程序中对应于Fortran的下标需要减一（比如，手册中提到，iparm(10)的参数，在C程序中，需要写为 iparm[9]）.