基于逆的多波前块ILU预处理方法

摘要

该发明公开了一种基于逆的多波前块ILU预处理方法，属于数值求解领域，目的是为了克服传统ILU预处理求解大型稀疏线性方程组时效率低下的问题。本发明主要采用了基于逆的抛弃策略，在很大程度上减少抛弃元素所引起的数值不稳定现象的发生，提高了数值稳定性；采用了不需存储更新阵的多波前法，避免了临时存储更新阵，大大减小了内存开销；采用了超块-自适应块不完全分解方法，使得计算性能、内存性能得到大幅提升。因此，本发明提出的基于逆的多波前块ILU预处理方法，可以大幅提升大型稀疏不对称线性方程组的求解效率，并且以很小的内存花销更加精确、快速地获得大型稀疏不对称线性方程组的解。

说明书

技术领域

本发明属于数值求解领域，具体涉及一种基于逆的多波前块ILU预处理方法，用于求解在科学和工程计算问题中产生的大型稀疏线性方程组。 

背景技术

在许多科学和工程领域，如核物理、高维微分方程数值解、结构和非结构问题的有限元分析、计算流体力学、石油地震数据处理、电力系统的优化设计及数值天气预报等，这些大型、超大型的科学与工程问题非常复杂。这时传统的实验方法、尝试方法不仅代价巨大，而且需要花费大量的时间；由于很难给出这些问题中各种物理量间的函数关系，无法通过解析方法得到这些复杂问题的精确解。相对于解析方法，数值方法并不通过具体数学函数关系求解，而是通过分布在若干点上的一系列具体数值来近似表征问题的解。数值方法虽只是近似方法，但只要所取点数足够多，也能以足够的精度逼近解析解，更重要的是它能够求解复杂的偏微分方程，因而能够在满足精度要求的情况下，解决解析方法不能解决的问题。特别是随着计算机技术的发展，诸如有限差分法、有限积分法、有限元法和边界元法等数值仿真方法被广泛应用于科学和工程问题，并成为最常用和最有效的方法之一。数值仿真方法的求解往往归结为大型稀疏线性方程组的求解，这也是此类方法中最消耗时间和内存的步骤。 

稀疏线性方程组的主要求解方法有直接法和迭代法两种，其中直接法通过直接对方程组进行一系列的运算而获得解；而迭代法通过逐次逼近来求解方程组，迭代直至达到预设精度为止。对于大型稀疏线性方程组，直接法存储花销和计算量巨大、计算成本高；而如果直接将迭代法用于求解，其收敛性和收敛速度得不到有效保证。因此，更为有效的方法是通过矩阵分裂预处理、不完全分解预处理以及稀疏近似逆预处理等方法首先对方程组进行预处理，再采用迭代法求解。其中，不完全分解预处理生成的矩阵由于能够对系数矩阵做较好的近似，能够大大的降低矩阵的条件数和迭代法的收敛步数、实现简单有效、适用范围广等优点得到了广泛应用。传统的不完全分解预处理有用于对称方程组的不完全cholesky分解预处理和用于非对称方程组的不完全LU(ILU)分解预处理等。 

传统的ILU分解预处理如ILUT分解(参见“Iterative Methods for Sparse Linear Systems”，Society for Industrial and Applied Mathematics，2003，Saad)在处理不到十万阶的小型稀疏线性方程组时具有很好的计算性能，而对于几十万甚至上百万阶的大型稀疏线性方程组，其计算性能并不让人满意。 

发明内容

本发明的目的是为了克服传统ILU预处理求解大型稀疏线性方程组时效率低下的问题，提出了基于逆的多波前块ILU预处理方法，利用该方法能够以更小的内存代价更加精确、高效地求出大型稀疏线性方程组的解。 

为了实现上述目的，本发明的技术方案是：基于逆的多波前块ILU预处理方法，包括以下步骤(除非特殊说明，本文中L和U分别表示下三角和上三角LU分解因子矩阵；P表示排序产生的排列矩阵)： 

步骤一：对原矩阵A执行重排序和符号分解，将矩阵重新组织成一系列的稠密矩阵， 

减少稀疏矩阵分解过程中产生的填入元数量，增加稠密操作。 

i).进行以减少填入元为目的的重排序； 

ii).执行符号分解，并建立消去树； 

iii).根据消去树后序遍历顺序进行重排序； 

iv).再次执行符号分解，建立消去树； 

最后形成重排序后的矩阵(其中P₀是重排序得到的排列矩阵，T表示转置)。 

步骤二：计算对角scaling矩阵D_r和D_c，得到矩阵增强数值稳定性。 

步骤三：将矩阵划分为若干超节点，以便利用稠密矩阵操作提升计算性能。 

步骤四：将各超节点划分为若干块矩阵。 

步骤五：采用基于超块-自适应块不完全分解方法和基于逆抛弃策略的不需存储更新阵的多波前法对各超节点进行块ILU分解(将第i个超节点F_i表示为块矩阵形式 $F_i=\left(\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right)$)。 

i).集成完全相加部分的值信息； 

ii).集成完全相加部分来自子孙超节点的更新； 

iii).执行完全相加部分块F₁₁的ILU分解其中是部分主元法产生的排列矩阵，是延迟主元产生的排列矩阵； 

iv).按照新的主元顺序调整已分解部分L₀,L₁,…L_i-1和U₀,U₁,…U_i-1，分别得到 和(其中下标0,1,…i-1表示第i个超节点之前的所有已分解的超节点)； 

v).按照新的主元顺序集成超节点的部分相加部分的值信息，得到调整主元顺序后的部分相加部分和

vi).计算并集成部分相加部分来自子孙节点的更新； 

vii).分别求解部分相加部分矩阵方程

获得F_i的块ILU分解所有超节点分解完毕后，可获得ILU分解 $P_r{P}_c\tilde{A}{P}_c^T=P_r{P}_c{D}_r{P}_{0}A{P}_0^T{D}_c{P}_c^T\approx \mathrm{LU}.$

步骤六：执行向前和向后回代，求解矩阵方程

i).执行向前回代，解(P_rP_cD_rP₀)^-1Ly＝b得

ii).执行向后回代，解Uz＝y得

iii).求解向量$x=P_0^T{D}_c{P}_c^{T}z.$

本发明的有益效果：利用本发明提出的基于逆的多波前块ILU预处理方法，可以大幅提升大型稀疏不对称线性方程组的求解效率，可以以很小的内存花销更加精确、快速地获得大型稀疏不对称线性方程组的解。本发明相对于传统的ILU预处理，在处理维数在十万阶以上的矩阵时，可以以2％～25％的内存优势获得10～18倍的速度提升。这是因为本发明提出的基于逆的多波前块ILU预处理方法有以下特点：1)采用了scaling技术，避免了一些数值难点，增强了数值稳定性并提高了计算精度；2)采用基于逆的抛弃策略，在很大程度上减少抛弃元素所引起的数值不稳定现象的发生，提高了数值稳定性；3)采用部分主元法、延迟主元法和对角扰动技术，进一步提高了鲁棒性和稳定性，使得预处理方法更加稳定可靠；4)采用不需存储更新阵的多波前法，避免了临时存储更新阵，大大减小了内存开销；5)采用超块-自适应块不完全分解方法，使得计算性能、内存性能得到大幅提升；6)所有块相关计算操作均采用BLAS3(Basic Linear Algebra Subprograms Level3，第三级基础线性代数子程序库)或LAPACK(Linear Algebra PACKage，线性代数函数库)中的稠密矩阵运算实现，大大提升了计算性能。 

附图说明

图1是本发明的主流程图； 

图2是矩阵的超节点划分图； 

图3是超节点的内部分块图； 

图4是波前阵分解流程图； 

图5是ILU分解流程图。 

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。 

如图1所示，基于逆的多波前块ILU预处理方法，包括以下步骤： 

步骤一：执行矩阵重排序和符号分解，将矩阵重新组织成一系列的稠密矩阵，减少稀疏矩阵分解过程中产生的填入元数量，增加稠密操作。 

所述步骤一包括： 

i).进行以减少填入元为目的的重排序 

采用基于多级嵌套剖分法和最小自由度排序技术的稀疏矩阵排序软件包METIS对原稀疏矩阵A进行重排序，以减少稀疏矩阵分解过程中产生的填入元数量，从而减小分解过程中的存储需求和计算量。 

ii).执行符号分解，并建立消去树 

执行符号分解，确定重排序后矩阵中的填入元位置，建立引导数值分解的消去树。 

iii).根据消去树后序遍历顺序进行重排序 

首先对消去树进行后序遍历，然后根据后序遍历顺序对矩阵进行重新排序，将矩阵重新组织成一系列的稠密矩阵，以获得更大的超节点，进而增大算法粒度、增加稠密操作。 

iv).再次执行符号分解，建立消去树 

再次执行符号分解，确定重排序后矩阵中的填入元位置，建立消去树。 

原矩阵A经过两次重排序后可表示为其中P₀是重排序矩阵，T是矩阵转置操作符。步骤二：计算对角scaling矩阵D_r和D_c，增强数值稳定性 

计算对角scaling矩阵D_r和D_c，使得scaling后的矩阵各行、列的范数近似为1。以避免大数除小数、值相近的数相减等数值难点，增强数值稳定性，加快分解速度，提高计算精度。 

步骤三：将矩阵划分为若干超节点，以便利用稠密矩阵操作提升计算性能 

采用超节点法将原矩阵划分为一系列超节点。超节点法将具有相同非零结构的连续列或行视作一个超节点，是多波前法中的一种常用的改进技术，以便利用稠密矩阵操作，减少了无效的间接寻址，增加了数据重用，可大幅提升算法的计算性能。由于系数矩阵具有对称的非零结构，划分出的超节点也是结构对称的。图2显示了系数矩阵的一种块划分结构，图中被划分为S₁,S₂,S₃,S₄四个超节点，超节点又由对角块矩阵(完全相加部分)和两 个在结构上对称的非对角块矩阵(部分相加部分)构成，其中最后一个超节点S4只包含对角块矩阵。 

为获得更大的超节点，进一步提高计算性能，本步骤中采用了以下两种松弛超节点技术： 

i).松弛超节点法：通过放松超节点非零结构“相同”的划分条件，允许在非零结构有一定比例的不同，通过显式添加零元获得更大的超节点； 

ii).超节点融合技术：将较小的连续叶子节点合并为一个超节点，限制最小超节点的维数。 

步骤四：将各超节点划分为若干块矩阵 

由于自然存在的块结构已经在系数矩阵的重排序过程中受到破坏，需要进行人为分块：将各超节点分别划分为若干维数固定的方块，在超节点边缘的块维数根据超节点的结构进行适当调整。超节点内部的分块情况如图3所示，a为设定方块维数，b和c是为了根据超节点完全相加部分和部分相加部分的边缘调整后产生的两种边长，图中超节点由八种不同尺寸的块构成。 

步骤五：采用基于超块-自适应块不完全分解方法和基于逆抛弃策略的不需存储更新阵的多波前法对各超节点进行块ILU分解 

不需存储更新阵的多波前法中，波前阵由超节点和其子孙超节点的更新阵通过矩阵扩充相加运算获得。设是来自第i个波前阵F_i的k个子孙节点更新阵，A_i是第i个超节点，则波前阵F_i满足如下关系 

$F_i=A_i+U_0^{\left(i\right)}+U_2^{\left(i\right)}+···+U_k^{\left(i\right)}---\left(1\right)$

其中运算符“+”代表矩阵扩充相加运算。由于超节点与其子孙节点具有不同的非零结构，需要通过矩阵扩充相加运算来获得波前阵。多波前法中波前阵的分解可表示为： 

$\left(\begin{array}{c}F=\left(\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& F_{22}\end{array}\right)=\left(\begin{array}{cc}{F}_{11}& F_{12}\\ F_{21}& 0\end{array}\right)+\left(\begin{array}{cc}0& 0\\ 0& F_{22}\end{array}\right)\\ =\left(\begin{array}{cc}{L}_{11}& 0\\ L_{21}& 0\end{array}\right)\left(\begin{array}{cc}{U}_{11}& U_{12}\\ 0& 0\end{array}\right)+\left(\begin{array}{cc}0& 0\\ 0& L_{21}{U}_{12}\end{array}\right)\end{array}\right)---\left(2\right)$

其中，L₁₁、L₂₁、U₁₁和U₁₂表示波前阵中相应块矩阵LU分解所得分解因子矩阵。传统的多波前法在分解完波前阵后，将计算当前波前阵对其父节点的更新阵，并将其保存至其父节点波前形成时才使用并释放，因此整个多波前数值分解过程中需要临时存储许多波前阵，这将消耗大量的存储资源。不需存储更新阵的多波前法中波前阵的分解由两个步骤组成：i)计算波前阵自身的分解：首先计算完全相加部分矩阵块ILU分解然后求解部分相加部 分矩阵方程(其中是部分主元法产生的排列矩阵，是延迟主元产生的排列矩阵)；ii)计算当前波前对其对祖先节点的更新L₂₁U₁₂，这个更新将被延迟到相应祖先节点波前阵分解时才计算并集成，因而避免了临时存储更新阵，大大减小了多波前法的内存需求。并因此避免了传统多波前法中在步骤一通过对消去树后序遍历顺序的优化来减少存储更新阵的操作，加速了步骤一的实现。 

由于完全相加部分主元顺序在分解过程中可能会发生变化，其余部分需要作出相应的调整。因此为减少主元顺序变化引起的重排序操作，不需存储更新阵的多波前法相应的主元顺序确定后才集成相应的之信息，例如：在完全相加部分分解完毕才按照主元顺序集成部分相加部分的值信息。此外，(2)式表明了部分相加部分两个块矩阵的的运算是无关的，可以单独进行，这里进一步将部分相加部分的分解分为上、下三角两个部分，波前阵的分解则变为了对完全相加部分块矩阵和部分相加部分块矩阵等相对较小的子矩阵分解，存储当前波前阵的临时存储需求进一步减小。 

在不完全分解过程中，抛弃元素可能导致对角元很小甚至为零，引起数值不稳定现象。基于逆的抛弃策略的采用将在很大程度上减少这种因素引起的数值不稳定现象的发生。对于LU分解A＝LU，预处理后的矩阵可表示为 

$\left(\begin{array}{c}{M}^{-1}A={\left(\tilde{L}\tilde{U}\right)}^{-1}A={\tilde{U}}^{-1}{\tilde{L}}^{-1}A\\ =(U^{-1}+Y)(L^{-1}+X)\mathrm{LU}\\ =I+\mathrm{YU}+U^{-1}\mathrm{XA}+\mathrm{YXA}\end{array}\right)---\left(3\right)$

其中和为相应的ILU分解的逆分解因子，逆分解因子的误差矩阵X和Y决定了预处理矩阵对原矩阵的近似程度。因此，基于逆的抛弃策略在门限抛弃策略的基础上，在抛弃条件中加入逆分解因子的范数，形成新的抛弃条件： 

$\left|l_{\mathrm{jk}}\right|{\left|\right|e_k^T{\tilde{L}}^{-1}\left|\right|}_{\infty }\le \tau ,\left|u_{\mathrm{km}}\right|{\left|\right|{\tilde{u}}^{-1}{e}_k\left|\right|}_{\infty }\le \tau ---\left(4\right)$

其中j＞k,m＞k，l_jk和u_km分别是分解因子L的第k列元素和U的第k行元素，单位向量e_k表示取矩阵第k行和第k列，τ为抛弃容限，T是矩阵转置操作符，分别是逆分解因子的第k行范数和的第k列范数。 

如图4所示，所述步骤五：包括： 

i).集成完全相加部分的值信息 

将超节点完全相加部分的值信息集成至波前阵。 

ii).集成完全相加部分来自子孙超节点的更新 

采用不需存储更新阵的多波前法，利用矩阵扩充相加运算将来自子孙节点的对完全相加部分的更新叠加至当前波前阵。 

iii).执行完全相加部分块ILU分解

首先估计分解因子的逆范数，然后确定抛弃条件； 

采用超块-自适应块不完全分解方法执行块ILU分解，这种方法包括三个方面：i)在抛弃方式上，它的基本抛弃单元是块，块中满足抛弃条件的元素都将被抛弃，若块中所有元素均被抛弃，将整个块抛弃；ii)在存储策略上，采用基于内存池的自适应存储策略，通过比较稀疏存储方式和稠密存储方式所需代价，自动选择存储代价小的存储方式；所有块存储空间和小片内存均采用内存池管理；iii)在计算架构上，采用超块方法，将与同一块或向量参与运算的多个连续存储的块拷贝到临时存储空间中形成一个超块，再利用稠密矩阵运算实现相关计算操作；iv)所有块的计算操作均采用BLAS3中的稠密矩阵运算实现，大大提高了预处理的计算性能。 

如图5所示，执行ILU分解时，首先采用部分主元法对完全相加部分进行分解，遭遇零主元时，采用延迟主元法将当前主元延迟消去，如此循环，直至分解成功。如果在所有主元顺序被尝试过后，依旧存在对角线元素过小或为零的或者分解失败，采用对角扰动技术，在对角线上添加扰动，以保证分解顺利完成。 

iv).按照主元顺序调整已分解部分L₀,L₁,…L_i-1和U₀,U₁,…U_i-1，分别得到和 

由于完全相加部分在分解时进行了行列交换，需要按照新的主元顺序对已分解部分作出相应调整。 

v).按照主元顺序集成超节点的部分相加部分的值信息 

按照新的主元顺序集成部分相加部分的值信息，得到调整主元顺序后的部分相加部分 和

vi).计算并集成部分相加部分来自子孙节点的更新 

利用矩阵扩充相加运算将来自子孙节点的对部分相加部分的更新叠加至当前波前阵。 

vii).分别求解部分相加部分矩阵方程和

至此，获得F_i的ILU分解所有超节点分解完毕后，整个数值分解阶段结 束，可获得ILU分解： 

$P_r{P}_c\tilde{A}{P}_c^T=P_r{P}_c{D}_r{P}_{0}A{P}_0^T{D}_c{P}_c^T\approx \mathrm{LU}---\left(5\right)$

步骤六：执行向前和向后回代，求解矩阵方程

用分解因子表示原系数矩阵有： 

$A\approx {\left(P_r{P}_c{D}_r{P}_0\right)}^{-1}\mathrm{LU}{\left(P_0^T{D}_c{P}_c^T\right)}^{-1}---\left(6\right)$

则线性方程组可表示为 

${\left(P_r{P}_c{D}_r{P}_0\right)}^{-1}\mathrm{LU}{\left(P_0^T{D}_c{P}_c^T\right)}^{-1}x=b---\left(7\right)$

i).执行向前回代，解(P_rP_cD_rP₀)^-1Ly＝b得

ii).执行向后回代，解Uz＝y得

iii).求解向量$x=P_0^T{D}_c{P}_c^{T}z$ 。