域分解算法作为求解器与预条件器的实现探究

在学习并行域分解（DD）概念并关注其实现时，“域分解作为求解器”与“域分解作为预条件器”是高频出现的术语。《并行计算机上偏微分方程的数值解法》（Bruaset 2006）中提到：

由于第1、2章介绍的施瓦茨方法（Schwarz methods）是应用于预条件化全局问题的不动点迭代，因此无法实现最快收敛速度，自然可以改用克雷洛夫方法（Krylov methods）。这就为将施瓦茨方法用作预条件器而非求解器提供了依据。

作为求解器，并行施瓦茨方法的逻辑十分清晰：只需在各子域上并行求解，然后交换边界数据即可。

但核心疑问在于：并行DD方法如何作为共轭梯度法（conjugate gradient method）的预条件器使用？预条件共轭梯度算法中，每一步应用预条件器时的全局同步是否会成为主要计算瓶颈？

尝试分析 #

核心思路是在各子域上进行局部求解，然后对这些解进行归约求和（reduce-sum）得到预条件器（此时可保存预条件器M_inverse），用于全局问题共轭梯度算法的每一步。但在更复杂的算法中，例如约束平衡域分解法（BDDC），共轭梯度算法的每一步都必须进行全局同步。此外，求解局部问题时还可为这些线性求解器添加预条件器，形成嵌套预条件算法（详见文末伪代码）。

《域分解方法导论：算法、理论与并行实现》（Dolean 2016）中限制加性施瓦茨（RAS）方法的定义如下：

以下是主进程控制下，将RAS作为预条件共轭梯度法（PCG）预条件器的伪代码：

def RAS(A_global, b_global, subdomain_indices):
    # 出自《域分解方法导论：算法、理论与并行实现》（Dolean 2016）Listing 3.2
    # 获取子域信息
    subdomain_id = get_partition_id() # 使用MPI
    A_subdomain = A_global[subdomain_indices[subdomain_id]]
    b_subdomain = b_global[subdomain_indices[subdomain_id]]

    # 局部求解也可添加预条件器
    jacobi_preconditioner_subdomain = diagonal(A_subdomain)
    x_subdomain = linear_solver(
        A_subdomain, b_subdomain, preconditioner=jacobi_preconditioner_subdomain)

    # 归约求和各子域解
    M_inverse = reduce_sum(x_subdomain) # MPI全局同步步骤
    return M_inverse


def preconditioned_conjugate_gradient(
    A_global, b_global, x0, subdomain_indices, n_iters):
    # PCG的收敛/迭代检查在主进程执行，预条件器为RAS
    # 出自《Scientific Computing An Introductory Survey》（Heath 2018）Algorithm 11.2
    # 初始残差
    rk = b_global - A_global @ x0

    # 通过RAS获取预条件器并保存，用于PCG每一步
    # 注意：并非总能保存M_inverse，每一步直接应用预条件器是更通用的实现！
    M_inverse = RAS(A_global, b_global, subdomain_indices)

    # 对残差应用预条件器
    sk = M_inverse @ rk

    # CG迭代
    for k in range(n_iters):
        alpha_k = (rk.T @ sk)/(sk.T @ A @ sk)
        
        # 更新x_{k+1}
        xk += alpha_k*sk

        # 更新r_{k+1}
        rkp1 = rk - alpha_k*A @ sk

        betakp1 = (rkp1.T @ M_inverse @ rkp1)/(rk.T @ M_inverse @ rk)

        # 应用预条件器得到s_{k+1}
        sk = M_inverse @ rkp1 + betakp1*sk 

        # 更新残差
        rk = rkp1

需要注意的是，实际应用中通常不会直接计算M_inverse，而是计算预条件器M对残差r的作用（即application/action），详见《Acceleration of a parallel BDDC solver by using graphics processing units on subdomains》（Sistek 2023）。《Numerical Linear Algebra with Julia》（Darve 2021）中的PCG方法伪代码更准确地展示了这一点：

def pcg(A, b, maxiter, preconditioner = BDDC):
    # 出自《Numerical Linear Algebra with Julia》（Darve 2021）第275、277、320页
    n = A.shape[0]
    x = zeros(n)
    r = b
    p = zeros(n)
    rho = 1.0
    for i in range(1, maxiter+1):
        # 应用预条件器：以BDDC为例，返回多个修正操作的和，即z = v1 + v2 + v3，各分量基于残差r计算
        z = preconditioner.apply(r)
        # 后续迭代步骤省略
        ...
        ...

    return x 

内容的提问来源于stack exchange，提问作者Jared