转载:Intel MKL 稀疏矩阵求解PARDISO 函数
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]).