将域分解算法实现为求解器与预处理器

Question

问题

在学习并行域分解 (DD) 概念以及随后对其实现感兴趣时，域分解作为求解器与预条件器这两个术语经常出现，请参阅下面的“并行计算机上偏微分方程的数值解”（布鲁塞特 2006）：

自从第 1 章和第 1 章介绍了 Schwarz 方法 2 表示应用于预处理全局问题的定点迭代，因此不提供可能的最快收敛，很自然地应用 Krylov 方法代替。这为使用 Schwarz 方法作为预处理器而不是求解器提供了理由。

作为求解器，至少并行 Schwarz 方法作为 DD 是相当清楚的，请参见this演讲中的下图，因为人们只是在每个域上并行求解，然后交换边界数据。 但是并行 DD 方法如何用作共轭梯度法的预处理器呢？作为在预条件共轭梯度算法中每一步应用预条件器的一部分的全局同步难道不会成为主要的计算瓶颈吗？

尝试

我认为这个想法是，你在每个子域上进行本地求解，然后对这些解决方案进行约简求和以获得作用于共轭梯度算法的每个步骤的预处理器（在这种情况下，可以保存预处理器

M_inverse

）全球问题。 然而，这似乎很奇怪，因为在更复杂的算法中，例如，通过约束平衡域分解（BDDC），必须在共轭梯度算法的每一步进行全局同步。最后，在解决局部问题时，您实际上可以对这些线性求解器使用预处理器，因此您最终会得到一种嵌套预处理算法（请参阅末尾的伪代码，可能会澄清这一点）。

限制加性 Schwarz (RAS) 方法的定义来自《领域简介》 分解方法算法、理论和并行实现”（Dolean 2016）：

RAS 的伪代码作为由主进程控制的预条件共轭梯度 (PCG) 方法的预条件器：

def RAS(A_global, b_global, subdomain_indices):
  # Listing 3.2 from "An Introduction to Domain Decomposition Methods Algorithms, 
  # Theory, and Parallel Implementation" (Dolean 2016)

  # Get subdomain info
  subdomain_id = get_partition_id() # using MPI
  A_subdomain = A_global[subdomain_indices[subdomain_id]]
  b_subdomain = b_global[subdomain_indices[subdomain_id]]

  # local solve can ALSO have a preconditioner
  jacobi_preconditioner_subdomain = diagonal(A_subdomain)
  x_subdomain = linear_solver(
    A_subdomain, b_subdomain, preconditioner=jacobi_preconditioner_subdomain)

  # add the solutions on each subdomain 
  M_inverse = reduce_sum(x_subdomain) # using MPI, this is a global sync step

  return M_inverse


def preconditioned_conjugate_gradient(
  A_global, b_global, x0, subdomain_indices, n_iters):
  # Convergence/iteration checks of pcg occurs on the master process, 
  # and its preconditioner is the RAS
  # Algorithm 11.2 from "Scientific Computing An Introductory Survey" (Heath 2018)

  # initial residual
  rk = b_global - A_global @ x0

  # get the preconditioner via RAS and save this so that it may 
  # be used at each step of the pcg algo, SAVING M_INVERSE
  # IS NOT ALWAYS POSSIBLE AND APPLICATION OF PRECONDITIONER AT EACH
  # STEP IS MORE GENERAL ALGO!
  M_inverse = RAS(A_global, b_global, subdomain_indices)

  # APPLICATION OF PRECONDITIONER ON RESIDUAL!
  sk = M_inverse @ rk

  # cg iterations
  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)

    # APPLY PRECONDITIONER!!! s_{k+1}
    sk = M_inverse @ rkp1 + betakp1*sk 

    # update rk
    rk = rkp1

值得注意的是，实际中

M_inverse

一般不是直接计算，而是计算预条件子M对残差r的

操作/应用

。下面以伪代码形式提供的“Numerical Linear Algebra with Julia”（Darve 2021）中的 PCG 方法更准确地说明了这一点。另请参阅 python

BDDC

实现和相应的 Python PCG。

def pcg(A, b, maxiter, preconditioner = BDDC):
  # PCG pseudo code with preconditioner application from
  # pg. 275, 277, and 320 from Darve 2021
  n = A.shape[0]
  x = zeros(n)
  r = b
  p = zeros(n)
  rho = 1.0
  for i in range(1, maxiter+1)
    # THE PRECONDITIONER! In BDDC case, this returns the sum
    # of different correction operations such that 
    # z = v1 + v2 + v3 where the components of the sum use
    # the residual r 
    z = preconditioner.apply(r)
    ...
    ...

  return x 

Answer 1

AMfeti

获得，我认为使用域分解方法作为预处理器与求解器本质上可以归结为这样一个事实：将域分割成不同的部分（称为子结构）后，您想要获得一个预处理器（即，结果）对于每个子域，预处理矩阵 s

 和残差 

M

（如 

r

）之间的矩阵向量乘法的 s = inv(M) @ r

，在该子域上使用迭代求解器，然后确保子域边界处的解不存在不以某种方式“冲突”。似乎在子域边界获得解决方案不“冲突”是域分解方法的不同之处（例如，FETI 使用拉格朗日乘子）。因此，本质上，域分解方法本身可以是求解器，或者域分解可以在每个子域上创建一个预处理器，以便与迭代求解器结合使用（例如，预处理共轭梯度）。我以图形方式总结了使用域分解方法作为每个子域的预处理器，然后单独的迭代求解器生成解决方案

u_{i}

，并随后解决边界冲突。解决边界冲突会在整个域上强制执行全局解决方案，并且此过程将持续下去，直到满足某些收敛标准。

最后一次验证预条件子在每个子域上都不同（在分布式设置中，每个子域可能存在于单独的计算节点上）这一事实得到了“C 和 Python 中偏微分方程数值解的 PETSc”的证实（Buehler 2021） ，我添加此内容只是为了在这个答案中显示我的主张的具体验证。下面的摘录分别有

bjacobi

表示块雅可比法和

asm

表示加法 Schwarz 方法：