[发明专利]一种高精度PSPIKE+求解器

权利要求书

1.一种高精度PSPIKE+求解器，其特征是该方法包括如下步骤：

(1)对线性系统矩阵进行重排和分割

高精度PSPIKE+求解器需要优化系统矩阵的对角占优性质以及系统矩阵的带宽，以提高求解器的求解精度和求解速度。其中，系统矩阵的对角占优性质会直接影响三角分解(LU分解)结果的精度。同时，系统矩阵的带宽越小，各个子系统之间的重叠区域越小，计算机的资源消耗越低。为实现矩阵对角增强及带宽降低，按照“先对角增强，后带宽降低”的顺序，分别对系统矩阵和右手边进行重排。为了对重排前后的线性系统做区分，假设重排前的线性系统为Fe＝h，重排后的线性系统为Ax＝b。重排后，根据带宽和分区数量，将线性系统分割为p个互相重叠的线性子系统，分割后第i个子系统通常表示为

其中，Aⁱ为第i个子系统的系统矩阵；xⁱ为第i个子系统的求解结果；bⁱ为第i个子系统的右手边；yⁱ_t和yⁱ_b为第i个子系统右手边中的顶部重叠区域和底部重叠区域的未知向量，用于表示第i-1个子系统与第i+1个子系统对第i个子系统的影响。

(2)计算耦合矩阵与逆矩阵

为提高求解器的求解速度，引入表示以及第i个与第i-1个子系统关系的耦合矩阵Alⁱ，第i个与第i+1个子系统关系的耦合矩阵Arⁱ。

其中，矩阵Aⁱ₂₂、Aⁱ₂₁和Aⁱ₂₃分别代表第i个子系统矩阵的对应部分。

对Aⁱ求逆，得到逆矩阵

高精度PSPIKE+求解器仅需要计算逆矩阵的一部分，包括矩阵Bⁱ₁₁、Bⁱ₁₃、Bⁱ₃₁和Bⁱ₃₃。这些矩阵可以利用简化的求逆过程直接获得，简化的求逆过程详见具体实施方式。

(3)构建平衡系统

线性系统Ax＝b中第i-1、i和i+1个互相重叠的子系统如下所示

第i-1个子系统底部τⁱ(τⁱ表示向量xⁱ₁的长度)行与第i个子系统顶部τⁱ行重叠，且第i+1个子系统顶部τⁱ⁺¹行与第i个子系统底部τⁱ⁺¹行重叠，因此第i-1和i+1个子系统的结果会受到第i个子系统结果的影响。为消除第i个子系统对第i-1和i+1个子系统的影响，需要引入向量y^i-1_b和yⁱ⁺¹_t

其中，向量xⁱ₁的长度与向量y^i-1_b相同，向量xⁱ₃的长度与向量yⁱ⁺¹_t相同，向量xⁱ₂的长度与向量y^i-1_b和yⁱ⁺¹_t不同。因为上述公式中包含向量长度与y^i-1_b和yⁱ⁺¹_t不同的xⁱ₂，所以难以直接利用上述公式构建平衡系统。为消去向量xⁱ₂，对不包含向量yⁱ_b和yⁱ_t的第i个子系统进行变形

将变形后的第i个子系统代入线性系统Ax＝b中

结合向量y^i-1_b和yⁱ⁺¹_t的定义，可以得到向量y^i-1_b和yⁱ⁺¹_t不包含xⁱ₂的计算公式

上述公式中除却待求解向量y^i-1_b和yⁱ⁺¹_t，还包含未知向量xⁱ₁和xⁱ₃。而在平衡系统中，除待求解向量以外，不能存在其他未知向量。为此，在第i个子系统的两侧同时乘上逆矩阵Bⁱ，提取其中顶部τⁱ行和底部τⁱ⁺¹行，得到xⁱ₁和xⁱ₃的表达式

为提高计算速度，引入向量btⁱ_t和btⁱ_b，计算公式为

将向量xⁱ₁和xⁱ₃的表达式代入到向量y^i-1_b和yⁱ⁺¹_t的表达式中相应位置，可以得到仅包含未知向量yⁱ_b和yⁱ_t的公式

根据上述公式，可构建用于求解形如My＝g的平衡系统，其中向量y包含子向量yⁱ_b和yⁱ_t。

(4)求解平衡系统

平衡系统的求解器采用BICG迭代求解方法，并在传统BICG法的基础上，将运算M×pv_j和M^T×pv^t_j分别转换为对向量pv_j和pv^t_j执行的不包含矩阵M和M^T的运算。由于矩阵M和M^T由耦合矩阵和逆矩阵组成，而耦合矩阵和逆矩阵的计算过程占用大量系统内存。因此，不使用矩阵M和M^T可以避免耦合矩阵和逆矩阵的计算，进而降低求解器对内存的占用。

为了降低内存的占用，对运算M×pv_j和M^T×pv^t_j进行替换，替换的运算过程包含如下步骤：矩阵M和M^T由多个子矩阵拼接而成，不计算各个子矩阵，分别将矩阵M和M^T按行划分为p个子矩阵(第i个子矩阵包含2×τⁱ行)，除去第i个子矩阵中全为零的部分(不包含任何矩阵的部分)，得到矩阵Mⁱ和MTⁱ；矩阵M和M^T的各个子矩阵均可用包含耦合矩阵和逆矩阵的表达式表达，根据矩阵M和M^T各个子矩阵的表达式，结合矩阵Mⁱ和MTⁱ在M和M^T中的位置，获取矩阵Mⁱ和MTⁱ的表达式；分别将向量pv_j和pv^t_j分割为p组子向量yⁱ_b和yⁱ_t，并分别按向量y^i-1_t、y^i-1_b、yⁱ_t、yⁱ_b、yⁱ⁺¹_t和yⁱ⁺¹_b以及y^i-1_b、yⁱ_t、yⁱ_b和yⁱ⁺¹_t的顺序组成列向量ymⁱ和ymtⁱ；矩阵Mⁱ和MTⁱ表达式的右侧分别乘上向量ymⁱ和ymtⁱ，并利用表达式，计算Mⁱ×ymⁱ和MTⁱ×ymtⁱ的结果，得到向量Cyⁱ_t和Cyⁱ_b以及Ctyⁱ_t和Ctyⁱ_b；按先Cyⁱ_t后Cyⁱ_b或者先Ctyⁱ_t后Ctyⁱ_b的顺序，合并各组向量Cyⁱ_t和Cyⁱ_b以及Ctyⁱ_t和Ctyⁱ_b，合并结果即为乘法运算M×pv_j和M^T×pv^t_j的结果。其中，Mⁱ和MTⁱ的表达式以及向量Cyⁱ_t、Cyⁱ_b、Ctyⁱ_t和Ctyⁱ_b的定义与矩阵M(M^T)的具体结构有关，详见具体实施方式。

上述运算过程中仅使用矩阵M和M^T的子矩阵Mⁱ和MTⁱ的表达式，并未直接使用由耦合矩阵和逆矩阵组成的M和M^T。但是，Mⁱ×ymⁱ和MTⁱ×ymtⁱ的表达式中依旧包含耦合矩阵和逆矩阵的部分子矩阵(Arⁱ、Alⁱ、Bⁱ₁₁、Bⁱ₁₃、Bⁱ₃₁和Bⁱ₃₃)右乘向量的运算。为避免矩阵Arⁱ、Alⁱ、Bⁱ₁₁、Bⁱ₁₃、Bⁱ₃₁和Bⁱ₃₃的使用，利用如下公式，对于Mⁱ×ymⁱ和MTⁱ×ymtⁱ中耦合矩阵和逆矩阵右乘向量的运算过程进行替代

其中，假设yitrpⁱ为运算过程中，耦合矩阵和逆矩阵右乘的向量。

采用上述M×pv_j和M^T×pv^t_j替换运算的平衡系统迭代求解方法包含如下步骤：利用预设的初始结果或上一次的迭代求解结果，分别计算向量r_j、r^t_j、pv_j和pv^t_j；利用M×pv_j和M^T×pv^t_j的替换运算，计算矩阵M和M^T分别右乘向量pv_j和pv^t_j的结果；利用M×pv_j的运算结果，结合向量r_j，计算步长参数α；以pv_j为梯度方向，以α为步长，计算更准确的迭代求解结果；利用M×pv_j和M^T×pv^t_j的运算结果，在向量r_j和r^t_j的基础上，结合步长参数α_j+1，计算向量r_j+1和r^t_j+1；重复上述步骤，直到向量r_j+1的范数小于收敛阈值，或迭代次数已达上限，输出迭代次数j_last以及迭代求解结果y_{j_last}。

(5)并行求解子系统

为计算第i个子系统中右手边，将平衡系统的求解结果y_{j_last}分割为i组向量yⁱ_t和yⁱ_b(长度分别为τ^i-1和τⁱ)，并将yⁱ_t和yⁱ_b代入到第i个子系统右手边表达式中相应位置处进行计算。在得到第i个子系统右手边后，对第i个子系统进行求解，以获取求解结果xⁱ。其中，各子系统的求解过程相互独立，便于使用并行计算方法进行加速。

(6)合并子系统求解结果

为获取线性系统Ax＝b的求解结果x，对子系统的求解结果xⁱ进行合并。对于子系统结果之间的重叠部分，取平均值以降低误差。在合并求解结果xⁱ后，按“先带宽降低，后对角增强”的顺序，对合并结果进行反向重排，得到反向重排后的求解结果e。

(7)外置迭代求解

利用直接求解过程(步骤(1)到(6))可以得到运算F^-1×h的结果。为提高结果的准确性，将直接求解过程代入到外置BICG迭代求解器中。在含预处理过程的BICG迭代求解器中需要计算F^-1×r_j+1和(F^T)^-1×r^t_j+1。其中，将向量r_j+1作为h，代入直接求解过程中，可得F^-1×r_j+1的结果。如果在步骤(1)中先对矩阵F执行转置运算，再进行矩阵重排，然后顺序执行步骤(2)到(6)，可得(F^T)^-1×r^t_j+1的结果。

外置BICG迭代求解器的主要流程为：利用BICG迭代方法的计算公式，结合预设的初始结果或上一次的迭代求解结果，计算向量r_j+1、r^t_j+1、pv_j+1和pv^t_j+1；利用直接求解器，分别计算F^-1×r_j+1和(F^T)^-1×r^t_j+1的结果z_j+1和z^t_j+1；利用向量z_j+1、r^t_j+1、pv_j+1和矩阵F，计算BICG算法中的步长参数α_j+1；以pv_j为梯度方向，α_j+1为步长，计算更为准确的迭代求解结果；重复以上步骤，直到向量r_j+1的范数小于收敛阈值或超过最大迭代次数，停止迭代，输出迭代次数j_last及迭代求解结果。

2.基于商业数学计算软件Matlab R2019b及PARDISO 6.0.0实现高精度PSPIKE+求解器。根据权利要求1所述的一种高精度PSPIKE+求解器的实现方法，其特征在于：引入耦合矩阵，用于表征子系统内部关系；结合耦合矩阵及部分逆矩阵，提出一种新的平衡系统构建方法，用于提高求解结果的准确性；根据平衡系统的构建方法，提出一种新的迭代求解方式，以降低内存占用，并节约耦合矩阵及部分逆矩阵计算时间。