Efficient Domain Decomposition for the Helmholtz Equation on GPUs. 1st Amit Rotem
arXiv:2606.21695v1 [cs.DC] 19 Jun 2026
Department of Mathematics Virginia Tech Blacksburg, USA [email protected]
Abstract—The Helmholtz equation governs wave propagation in acoustics, electromagnetics, and seismology, but its indefinite nature makes it difficult to solve with iterative methods. Domain decomposition methods are a natural fit for massively parallel architectures, yet mapping efficient Helmholtz solvers onto modern GPUs remains a challenge. We address both with two key contributions: (1) a block-level domain decomposition scheme, in which each subdomain is assigned to a single thread block and all solves run concurrently in a single kernel launch, and (2) WaveHoltz as the subdomain solver. WaveHoltz is a fixedpoint iteration that is uniquely well-suited to the GPU execution model due to its minimal memory footprint and no reduction operations. Together, these eliminate device-level synchronizations and replace global memory traffic with shared memory and register-level operations, keeping subdomain data largely resident in L1 and L2 cache. We explore two threading strategies: one degree of freedom per thread for small subdomains, and multiple degrees of freedom per thread for larger ones. Benchmarks of our CUDA based implementation on a NVIDIA A100 show that WaveHoltz achieves 2x–25x speedup over MINRES, with the advantage growing with subdomain size. Crucially, evaluating the subdomain solver in single rather than double precision yields an additional 2x–10x speedup–a benefit largely unattainable by MINRES due to loss of Krylov vector orthogonality under reduced precision. Index Terms—High Performance Computing, GPU, Waves, Linear Algebra, FEM.
I. I NTRODUCTION Due to its indefinite nature, the Helmholtz equation is difficult to solve with iterative methods [1], [2]. Krylov space methods, such as MINRES/GMRES, converge very slowly without a good preconditioner. Other iterative methods such as an off-the-shelf multigrid do not converge for this problem and must be nontrivially specialized to the Helmholtz equation. Algorithms based on factorization, such as multifrontal methods, or ILU and sweeping preconditioners can be very effective, but the inherent sequentialism and superlinear costs of these methods makes them less suitable for massively parallel GPU computing. Time domain methods such as controllability methods [3], [4] and WaveHoltz [5], [6], [7] are also effective. This material is based upon work supported by the National Science Foundation under Grant Numbers DMS-2345225 and DMS-2436319 and Virginia Tech. This work has been submitted to IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible.
Fig. 1. A partition of a rectangular domain Ω ⊂ R2 into several subdomains. The elements are outlined in black, and the subdomains are outlined in purple. Two representative subdomains Ωℓ and Ωr and their common interface Γℓ,r are highlighted. The solution visualized is the wave produced by an off-center point source in a Y-shaped wave guide with a 3× contrast in wavespeed relative to the surrounding media. The solution is computed with degree 3 basis functions using the method described here.
The method we present here is a non-overlapping domain decomposition method. Domain decomposition methods (DDMs) are reliably efficient for solving the Helmholtz equation and naturally parallelize to many cores. The Helmholtz equation is indefinite and DDMs that solve it differ slightly from related methods for positive definite problems, mainly in the choice of transmission conditions between subdomains (discussed further in Section III.) The first method extending DDMs to the Helmholtz equation was presented in [8]. Recent advancements towards making DDMs more efficient typically target the convergence rate of the iteration. One class of methods aim to improve the transmission condition, for example, by incorporating perfectly matched layers [9], [10], [11],
[12], [13] or approximately absorbing boundary conditions (ABCs) [13], [14], [15], [16], [17]. The latter of which can be further improved by optimizing the coefficients of the ABC for the given problem [15]. Another class of methods improve the convergence rate of the iteration by a two-level coarse-space correction also known as deflation. The deflation vectors are computed either as part of the Krylov iteration [18], [19], [20] or by choosing a suitable ansatz space such as planewaves or eigenvectors of the transmission operator, see for example [21], [22], [23]. The convergence rate of DDMs scales adversely with the number of subdomains, but coarsespace methods effectively overcome this difficulty [20], [21]. The best performing DDMs typically employ both improved transmission conditions and coarse-space correction. In the context of high performance computing, domain decomposition naturally scale to many core environments. A detailed account is given in [24]. Typically, one subdomain is assigned to one node and communication takes place between nodes that share an interface. Compared to a Krylov method (without a preconditioner) for the same problem distributed the same way, DDMs tend to be very efficient not only because they take far fewer iterations but also because the work between communication stages is much greater. To make the latter statement true, the work per CPU must be large. For example, on the order of one hundred thousand DOFs per core were used in the experiments conducted in [24]. On GPUs, this requirement is pronounced. When used effectively a single GPU can process massive amounts of floating point operations compared to a single CPU node and therefore must be saturated with much more work to outweigh communication costs. To that end, there are two immediately obvious strategies to make domain decomposition strategies efficient on GPUs: (A) larger subdomains, or (B) multiple subdomains on one GPU. Strategy (A) is merited by the fact that DDMs with a few large subdomains will converge faster than those with many small subdomains. This comes at the cost of solving a larger Helmholtz equation on each subdomain which, as already discussed, is difficult and the cost may not scale optimally with the problem size. Here we consider strategy (B) which in turn can be used as an efficient inner solver for strategy (A). To make strategy (B) effective, we propose taking moderately sized subdomains consisting of several hundred to a few thousand degrees of freedom (DOFs) and solving each subdomain problem using a thread block consisting of enough threads so that each thread is associated with a small fixed number of DOFs. We present two variants of this strategy. First, we assign one thread to one DOF and the subdomains are taken small enough (fewer than 1024 DOFs) to store the entire subdomain in shared memory. The resulting kernel has a high arithmetic intensity and nearly all data movement occurs between registers and L1 cache. In the second variant, we have one thread per T DOFs (2 or 4, say), and use a smaller shared memory buffer together with a large global buffer. This variant can solve problems on larger subdomains (fewer than 1024 × T
DOFs) at the cost of increased cache and register pressure. The alternative to these strategies is to distribute each subdomain over several thread blocks. For large subdomains, doing so may be practically unavoidable but comes at the cost of device level synchronizations. For such a configuration, iterative methods are less optimal. Solving the subdomain problem may require hundreds of kernel launches compared to just one with the strategy we propose. When split over many kernels, the arithmetic intensity decreases and data movement dominates. If larger subdomains are ultimately necessary, then it may be more advantageous to employ different strategies such as direct factorization methods, or using our proposed method hierarchically, but we do not explore this here. Our approach is especially effective in two dimensions where subdomains with a few hundred degrees of freedom are large relative to the global problem, so the global system converges quickly. In three dimensions, these subdomains are notably smaller relative to the global problem, so other global acceleration strategies (such as coarse-space correction) are needed to make the solver effective. The paper is organized as follows. In Sections II–IV we introduce the Helmholtz equation, the non-overlapping domain decomposition method, and its spectral element discretization. In Section V we describe the WaveHoltz iteration for the subdomain solve, and in Section VI we combine these components into a GPU algorithm and propose two threading strategies. Finally, Section VII benchmarks the proposed method n two and three dimensions and compares WaveHoltz with MINRES as subdomain solvers. II. BACKGROUND Consider the Helmholtz equation in a bounded and simply connected domain Ω ⊂ Rd (d = 2, 3): −∆u − k 2 u = f,
x ∈ Ω,
(1a)
∂n u − iku = λ,
x ∈ ∂Ω.
(1b)
Here k = ω/c(x) where ω is the time-harmonic frequency and c(x) > 0 is the wave speed. The function f = f (x) is a compactly supported forcing term (typically real-valued) and λ = λ(x) is complex-valued. When λ = 0, the impedance boundary condition approximates the Sommerfeld radiation condition d−1
lim r 2 (∂r u − iku) = 0,
r→∞
r = |x|.
This radiation condition enforces the intuitive idea that all waves propagate out from the source f and that no waves are propagating from the far-field. More accurate approximations of the radiation conditions such as perfectly matched layers, or high order absorbing boundary conditions (see [13], [25], [26]) can meaningfully improve the quality of the solution, but here we limit our consideration to the impedance condition (1b). Nevertheless, it is important to note that such boundary conditions are particularly relevant to domain decomposition methods which discuss we further in Section III.
III. D OMAIN D ECOMPOSITION In this section, we broadly introduce the parallel nonoverlapping domain decomposition methods for the Helmholtz equation. Define LSnon-overlapping subdomains Ωℓ for ℓ = L 1, . . . , L such that ℓ=1 Ωℓ = Ω and Ωℓ ∩ Ωr = Γℓ,r = Γr,ℓ where Γℓ,r is a d − 1 dimensional interface such as in Figure 1. We define uℓ as the restriction of u to subdomain Ωℓ . The domain decomposition method is iterative, and at iteration n + 1 we update the solution un+1 by solving a Helmholtz ℓ equation on Ωℓ . Namely, −∆un+1 − k 2 un+1 = f, ℓ ℓ
IV. S PECTRAL E LEMENT D ISCRETIZATION We now present the weak formulation of (1). We assume that the domain is partitioned into convex quadrilateral or hexahedral elements which each element belonging to only one subdomain. We assume quadrilateral or hexahedral elements to exploit the tensor product structure of the basis functions. We seek u ∈ H 1 (Ω) such that for all ϕ ∈ H 1 (Ω)
x ∈ Ωℓ
Bℓ,r un+1 = Bℓ,r unr , ℓ (∂n − ik)un+1 = 0, ℓ
(∇u, ∇ϕ) − (k 2 u, ϕ) − ⟨ku, ϕ⟩ = (f, ϕ) + ⟨λ, ϕ⟩.
x ∈ Γℓ,r \ ∂Ω,
Enforcing the interface condition requires evaluating the normal derivatives of the solution at the boundary. This proves to be problematic for element based methods for which the derivatives are not necessarily continuous. To resolve this issue, we introduce the dual variables λℓ,r for each pair of neighboring subdomains (Ωℓ , Ωr ) such that on Γℓ,r
We then solve the subdomain problems with the interface conditions: x ∈ Γℓ,r .
It can be verified that from their definition (see [27]) that the iterates λnℓ,r satisfy: ℓ,r
a(x)b(x) dx,
(a, b) =
⟨a, b⟩ =
a(x)b(x) ds. ∂Ω
Ω
After discretizations by the spectral element method [28] we are left with: (A − iωH − ω 2 M )u = b.
(4)
Here the matrices A, M, H ∈ CN ×N and the vector b ∈ CN are defined as follows: Aij = (∇ϕj , ∇ϕi ),
Mij = (c−2 ϕj , ϕi ),
Hij = ⟨c−1 ϕj , ϕi ⟩,
bj = (f, ϕj ) + ⟨λ, ϕj ⟩,
where {ϕi }N i=1 is the basis for the spectral element method. It will prove useful in Section V to approximate the integrals with the GLL quadrature rule on which the nodal basis is defined resulting in M and H being diagonal (mass-lumped). We can also formulate (4) as a real block system for the real and imaginary parts of u: A − ω2 M ωH ℜ{u} ℜ{b} = . (5) ωH ω2 M − A ℑ{u} −ℑ{b} From the definitions of A, M, and H it is clear that the block system in (5) is symmetric, nevertheless it is indefinite making it difficult for iterative methods to converge. We also note that the matrices A, M, and H or the block system (5) need not be explicitly formed. The linear system can be solved iteratively with a Krylov method with matrix-vector products computed by evaluating integrals on the fly via quadrature. For more details, we refer to [29] and [30] wherein finite element operators are computed matrix-free on GPUs. V. S OLVING THE S UBDOMAIN P ROBLEM : WAVE H OLTZ
λr,ℓ = (∂n − ik)uℓ .
n n+1 λn+1 . ℓ,r = −λr,ℓ + 2ikur Γ
Z
Z
Bℓ,r = ∂n − ik.
(∂n − ik)un+1 = λnℓ,r , ℓ
(3)
Here,
x ∈ ∂Ω.
The transmission condition in terms of the operators Bℓ,r on Γℓ,r characterize the particular domain decomposition method and the choice of operator here critically affects the convergence rate of the iteration. In the classical parallel Schwarz method Bℓ,r = Id, but this choice generally leads to a non-convergent iteration for the Helmholtz equation and the subdomain problems may not be well posed [27]. The theoretically optimal choice is to take Bℓ,r = ∂n + DtN where DtN is the Dirichlet-to-Neumann (also known as the Poincaré-Steklov) operator. For a partition of Ω into L vertical slices, say, the iteration with the DtN transmission condition converges in exactly L iterations. However, the DtN is a non-local operator and cannot be discretized efficiently. There are many ways to approximate the DtN by local operators namely by ABCs just as discussed in Section II for the Sommerfeld radiation condition, see [9], [11], [13], [15]. A good choice of boundary conditions can dramatically improve the convergence rate [12]. In this paper, we are concerned moreso with the efficient implementation of these methods on GPUs, so for simplicity we take the so called zero-th order transmission condition
λℓ,r = (∂n − ik)ur ,
Note that rather than iterating on uℓ directly, we can equivalently iterate on λℓ,r instead. This fixed point iteration can be accelerated with a Krylov space method and the resulting method is equivalent to preconditioning the global system [27].
(2)
On each subdomain Ωℓ , (1) is discretized with the spectral element method as described in Section IV which is equivalent to solving (4) or (5). Among the fastest serial methods for solving this subdomain problem is to pre-compute a sparse factorization of (4) associated with each subdomain and perform the sparse triangular solves at every iteration. This method works very well in message passing parallelism because the work per node is large relative to the communication cost [24]. While there are many efforts to accelerate sparse triangular
solves on GPUs, for example [31] and packages like CuDSS [32], the process is inherently sequential. To take advantage of the highly parallel nature of the GPU, we propose that the linear system should be solved iteratively (hopefully in a small number of iterations). In Section VI we propose solving each subdomain problem on a single thread block consisting of several hundred threads, so the iterative method we use should be highly parallelizable with a small memory footprint as the storage will occupy registers and shared memory which are highly limiting resources on the GPU. The Krylov space methods MINRES and BiCGSTAB meet our requirements [33] and in Section VII we show performance results using MINRES as it is typically more reliable than BiCGSTAB and takes advantage of the symmetry of the problem at roughly the same cost. Here, we proceed with the WaveHoltz iteration [5] which has a smaller memory footprint than both MINRES or BiCGSTAB, requiring only one additional vector to carry out the iteration where MINRES requires five vectors. A second benefit of the WaveHoltz method is that it does not require any reduction operations (inner products) greatly reducing the number of thread synchronization steps. The WaveHoltz iteration solves the Helmholtz equation by filtering the solution of the wave equation in time. Starting from (4), let w = ue−iωt , then w solves: M ẅ + H ẇ + Aw = be−iωt , with initial conditions w(0) = u and ẇ(0) = −iωu. In fact, it is sufficient to solve for just the real part of w, say p = ℜ{w}, then M p̈ + H ṗ + Ap = ℜ{be−iωt }. (6) The WaveHoltz iteration solves for u by finding initial conditions to p so that p is periodic with frequency ω. We start with an initial guess u and perform the fixed point iteration (m+1)
u
(m)
= Πu
Nt −1 1 2 X 1 (κν pν − iω κν+ 12 ṗν+ 2 ). (7) = Nt ν=0
Here Nt is the number of time steps, κν = cos 2πν Nt − 1 1 2πν 2 π ν 4 + 4 tan Nt is the filter, and p = p( Nt ) where p(t) solves (6) with initial condition p(0) = ℜ{u(m) } and ṗ(0) = ℜ{−iωu(m) }. The filter κν selects for the ω-frequency component of p(t). The wave equation is evolved with the following modified leapfrog scheme. For the j-th component, h i ν− 1 ν+ 1 ṗj 2 = αj ṗj 2 + βj −(Apν )j + ℜ{bj e−2πiν/Nt } , ν+ 21
pν+1 = pνj + σ ṗj j Here αj = tan(ω∆t/2) ω
Mjj −θHjj Mjj +θHjj
and βj = sin(ω∆t/2) . ω/2
.
σ Mjj +θHjj
where θ =
and σ = Observe that θ ≈ ∆t/2 and σ ≈ ∆t are chosen so that if p0 = ℜ{u} and 1 ṗ 2 = ℜ{−iωue−iω∆t/2 } then pν = ℜ{ue−iνω∆t } and 1 1 ṗν+ 2 = ℜ{−iωue−i(ν+ 2 )ω∆t }. That is, this scheme is designed to exactly evolve e−iωt for fixed ω. This property, together with the time-filter κν guarantee that the WaveHoltz iteration exactly converges to u. Note also that this scheme
relies on M and H being diagonal so that the scheme is effectively explicit. This scheme has a time-stepping restriction ∆t ≤ p2 2h cmax where h is the smallest characteristic element length, p is the degree of the basis, and cmax is the fastest wavespeed in Ω. In [6], it was demonstrated that WaveHoltz can be interpreted as a preconditioned iteration for (4). For large problems the WaveHoltz iteration can be further accelerated with a Krylov space method and preconditioned with deflation making it robust and efficient solver on its own; see [7] where also implicit time-stepping is used. However, to reduce memory and syncronization cost, we do not accelerate the iteration in this work. We find that very few (five to twenty) iterations are sufficient to converge to machine precision for the sizes of subdomains in our method. One might also consider solving (4) on each subdomain to a low tolerance initially, then tightening the tolerance as (2) approaches convergence to reduce total runtime. We do not investigate this strategy here, but report that significant progress towards convergence can be made with as few as two WaveHoltz iterations. S UMMARY OF THE DDM C YCLE One step of the iteration (2) on one subdomain is summarized in Algorithm 1. It is clear from the pseudocode that to solve the subdomain problem, we need only store a complex vector u (or equivalently two real vectors) and two real vectors p and ṗ The operations on the vector v can be performed inplace on u. Similarly, Ap need not actually be stored, the result of the matrix-vector multiplication can be added directly to ṗ. In comparison, an efficient MINRES implementation will require five complex vectors and two inner products per iteration. We could monitor the convergence of WaveHoltz by checking the residual of the fixed-point iteration, but it is worth noting that far fewer inner products are needed by WaveHoltz compared to MINRES because WaveHoltz requires far fewer iterations at the cost of more work per iteration. VI. GPU A LGORITHMS In this section, we summarize the design of the domain decomposition method for the GPU. The ultimate goal is to exploit as much parallelism as possible in the domain decomposition method. At each outer iteration, we apply the DDM cycle of Algorithm 1. Since we can apply the DDM cycle to all subdomains concurrently, we identify this as the highest level parallelism we can exploit. Indeed this is the key feature that makes domain decomposition methods scalable [24], [27]. By solving these linear systems with an iterative method as described in Section V we expose much finer grained parallelism through the computation of vector operations and in the evaluation of SEM operators. Synchronizations are ultimately unavoidable in the evaluation of the SEM operators, namely, the evaluation of Ap. There are several levels of synchronizations that can take place on a GPU; warp-level, block-level, and device-level synchronizations all come with their own costs. Warp-level synchronizations are essentially free by nature of how the GPU
Input: λn Output: λn+1 Assemble b from f and λnℓ,r for each neighboring Ωr ; Initialize u(0) ← 0; while not converged ; // WaveHoltz do Initialize p ← ℜ{u(m) }; Initialize ṗ ← ℜ{−iωu(m) e−πi/Nt }; κ1 Initialize v ← κ0 p + 2 ṗ; iω for ν ← 1 to Nt − 1 ; // time stepping do foreach j do pj ← pj + σ ṗj ; end Compute Ap ; // Matrix-free foreach j do ṗj ← αj ṗj + βj −(Ap)j + ℜ{bj e−2πiν/Nt } ; end κν+ 12 ṗ; v ← v + κν p + iω end u(m+1) ← v; end n λn+1 r,ℓ ← −λℓ,r + 2ik u Γℓ,r ; Algorithm 1: (DDM Cycle) Domain decomposition subroutine on subdomain Ωℓ .
schedules instructions to threads. Block-level synchronizations are slower, but far slower are device-level synchronizations. Looking at the data dependence in Algorithm 1, we observe that synchronization must occur between WaveHoltz iterations, between time-stepping iterations, and during the computation of Ap. If the DDM cycle is implemented at the device level, that is, a single domain is treated by more threads than the maximum number of threads per thread block, then the number of independent kernel calls must scale like O(Nt ) per WaveHoltz iterations (for MINRES or a similar iterative solver, fewer kernel calls are needed per iteration, but many more iterations are needed). The arithmetic intensity per timestep is low and each kernel call must read and write to p and ṗ which necessarily occupy global memory. If each subdomain can be treated by a single thread block then this data movement is redundant. Instead we advocate that the entire DDM cycle be executed by a single thread block in a single kernel call. Then the number of reads and writes to global memory are substantially reduced and replaced with reads and writes to shared memory. To emphasize this point, note that NVIDIA’s A100 GPU used in our benchmarks in Section VII has a latency of approximately 300 cycles for global memory, 200 cycles for L2, and 20-30 cycles for L1/Shared memory as reported in [34]. Moreover since all compute operations are fused into one kernel the compute cycles dominate the runtime so the GPU is being effectively used.
To perform the DDM cycle on a thread block, we explore two threading strategies. First we assign one degree of freedom (DOF) to each thread. Second, we assign T DOFs to each thread with T small. In our examples we take T = 2 and 4. We now describe these strategies. A. One DOF per Thread On every NVIDIA GPU since 2010 only 1024 threads can be grouped into one thread block implying that subdomains of at most 1024 DOFs can be considered for our proposed strategy. To evaluate Ap matrix-free, it is more practical to describe the method in terms of elements rather than individual DOFs. To that end, assume a subdomain has K elements of degree P and assign to each thread block K(P + 1)d threads. While K(P + 1)d exceeds the number of DOFs, the redundancy in threads is worthwhile for computing Ap matrixfree as it is by far the most expensive operation. Note that the number of threads needed in three dimensions per element is substantially greater than in two dimensions limiting the number of elements that can be assigned to each subdomain. It is not necessary that each subdomain have exactly K elements. We assume that all subdomains have at most K elements and excess threads go unused1 . At the start of Algorithm 1, all threads read the right hand side f from global memory, and the threads associated with boundary DOFs additionally read the transmitted values λ. Each thread stores a single entry bj assembled from f and λ independently in registers. The threads store single entries of p, ṗ, and u and their corresponding entry of α and β in registers. Note that all steps in Algorithm 1 besides the matrix-vector product Ap are entry-wise operations that can be carried out by the threads concurrently without any branching, synchronizations, or data movement outside of registers. Computing Ap requires synchronization and the sharing of memory between threads. The shared memory consists of three buffers, the 1D differentiation matrix D of size (P + 1)2 , a buffer pShared of length K(P + 1)d where the solution is stored, a buffer gradShared of length dK(P + 1)d where the gradient is stored. Shared memory usage can be further reduced at the cost of one additional synchronization per evaluation of Ap if pShared aliases gradShared. When evaluating Ap, each thread maps its local entry of p to its corresponding spot in pShared and the threads synchronize. Each thread computes the gradient of p at one quadrature point, which consists of d inner products with the corresponding rows of D. The threads write to gradShared then synchronize. Each thread accumulates the integral (∇p, ∇ϕ) which consists of another d inner products with rows of DT . The result is accumulated in pShared atomically and the threads synchronize. The atomic operations are necessary since DOFs on element faces accumulate contributions from more than one element and therefore more than one thread. This does not incur any substantial serialization 1 To avoid warp divergence, the DOF vector is extended by zeros and the excess threads accumulate zeros.
because only a few threads write to the same location; for the most part, threads write independently. After the WaveHoltz iteration converges, the threads associated to boundary DOFs compute λn+1 and write to global memory. This summarizes the one-DOF-per-thread approach. B. T DOFs per Thread
VII. B ENCHMARKS A. Efficiency in Two Dimensions To evaluate the performance of the parallel strategies described, we measure the runtime of the DDM cycle with five WaveHoltz iterations applied to a random input vector. We compare the runtime relative to the time it takes to evaluate the Helmholtz operator A−iωH −ω 2 M on the whole domain applied to the same random input vector. We consider a constant coefficient problem in [−1, 1]2 with uniform elements.
2. Map a chunk of elements to shared memory buffer.
3. Threads compute chunks of 𝐴𝐩 in shared memory and store partial result in registers Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Shared Memory
Global Memory
Now assume each subdomain has no more than T × K elements and assign to each thread block K(P + 1)d threads. We assume that T K(P +1)d > 1024 or that a buffer of length dT K(P + 1)d does not fit in shared memory. For the most part, the procedure is the same for the one DOF per thread case, but each thread now stores T entries of b, p, ṗ, and u and the entry-wise operations in the DDM cycle are serialized by a factor T compared to the one DOF per thread approach. The key difference between the two approaches is in the computation of Ap. As before we have the shared memory buffers for D, pShared consisting of K(P + 1)d entries, and gradShared consisting of dK(P + 1)d entries. Note that this requires T times fewer entries than above relative to the number of DOFs. Additionally, a global buffer pGlobal of size L × T K(P + 1)d is needed (roughly the same size as the solution vector) and each subdomain operates on a slice of length T K(P + 1)d of this buffer. When evaluating Ap, each thread maps all T values of p to pGlobal and the threads synchronize. The threads now work to evaluate Ap in blocks of K elements concurrently. For each block of elements, the threads copy K elements from pGlobal to pShared. The threads compute the entries of Ap associated with this block of elements like in Section VI-A but store the element contributions in registers rather than accumulate in pShared. The threads repeat this process T times to cover all T K elements resulting in T entries of Ap in registers. Then the threads accumulate the entries of Ap in pGlobal atomically and synchronize. Finally, the threads collect the accumulated values of Ap back from pGlobal into registers. This procedure is visualized in Figure 2. In Section VI, we emphasized that data movement to global memory can take as long as 300 cycles which suggests that in the worst case, the T DOFs per thread strategy may be substantially more expensive than the one DOF per thread strategy. However, due to the frequent and predictable access pattern of the memory access to a fairly small slice of global memory (no more than 1024T values), the cache can effectively amortize the cost of the memory movement which we show in Section VII.
1. Map register values to global work buffer.
…
6. Map results from global back to registers.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
Reg.
4. Synchronize and repeat for each chunk of 𝐴𝐩
5. Accumulate chunks of 𝐴𝐩 in global work buffer.
Fig. 2. The T DOFs per Thread algorithm with six threads and four DOFs/thread. Arrows indicate memory movement. When T = 1 (Section VI-A), the memory traffic through global is bypassed, and only the shared memory buffer is used.
Since the DDM is a preconditioner, we can gain significant speedups simply by evaluating the preconditioner in a lower precision, for example as considered in [35], [36], [37]. In our numerical experiments we observed roughly 1.5× to 3.5× speed up to be gained by evaluating the DDM cycle in single precision on an A100. Some GPU architectures such as Ada Lovelace emulate 64 bit floating point numbers and therefore achieve much lower 64 bit throughput relative to 32 bit throughput. We measured a 5× to 10× times speedup on an RTX 4060Ti. The following experiments are performed in single precision on an A100. We evaluate the operators for a sequence of meshes from 32 × 32 elements to 1024 × 1024 elements of polynomial degree P ∈ {1, 2, 3, 4, 7}. We vary the thread block sizes B ∈ {256, 512, 1024} and the DOFs per thread T ∈ {1, 2, 4}. For T = 1 we employ the strategy discussed in VI-A and for T = 2 and 4 we employ the strategy discussed in VI-B. For each combination of P, B, and T we set the number of elements per subdomain as large as possible. For example, for P = 3, B = 1024, and T = 2 we take ⌊ (PBT +1)2 ⌋ = 128 elements in each subdomain, e.g. as blocks of 16×8 elements. In Figure 3 we plot the relative time for three kernels. Observe that the DDM is substantially more expensive for high order than lower order discretizations. We attribute this scaling to the time stepping restriction ∆t ≲ 2h/p2 . The cost of the WaveHoltz solver scales directly with 1/∆t, so we expect
B = 256, T = 1
50
B = 1024, T = 1
50
30
DD Time / Operator Time
40
DD Time / Operator Time
40
DD Time / Operator Time
40
30
20
30
20
10
20
10
1×
0
104
105
#DOFs
106
10
1×
0
107
B = 1024, T = 4
50
104
P=1
105
#DOFs
106
Degree P=3
P=2
1×
0
107
P=4
104
105
#DOFs
107
106
P=7
Fig. 3. The time to evaluate the DDM cycle (DD time) relative to the time it takes to evaluate A − iωH − ω 2 M on the entire domain (Operator time) for a sequence of problems of increasing size in two dimensions. P=3
P=1
Throughput [MDOF/s]
Throughput [MDOF/s]
Throughput [MDOF/s]
500
400
400
400
300
300
300
200
200
200
100
100 0
P=7
500
500
104
105
#DOFs
B = 256, T = 1 B = 256, T = 2
106
0
B = 256, T = 4 B = 512, T = 1
100
106
105
#DOFs
Kernel B = 512, T = 2 B = 512, T = 4
107
0
B = 1024, T = 1 B = 1024, T = 2
105
106
#DOFs
107
B = 1024, T = 4
Fig. 4. The throughput (millions of DOFs per second) of the various kernels for fixed degree and increasing problem size.
the DDM to be more expensive for higher order discretizations2 . As the problem size increases, the cost ratio tends to depend largely on degree and less so on kernel configuration. Nevertheless, kernels operating on larger subdomains appear to be less efficient for smaller problems. For example, the B = 1024, T = 4 kernel is more efficient for larger problems than smaller problem for all degrees. Notably, the DDM cycle is very efficient for P = 1 and can be evaluated only 2× to 4× the cost of the global operator. In Figure 4, we plot the throughput in millions of DOFs per second for p = 1, 3, and 7 for all combinations of B and T . For p = 1 and p = 3 the kernels operating on the smallest subdomains are most efficient. In contrast, the most efficient kernels for p = 7 are B = 256, T = 2 and B = 1024, T = 2 and T = 4. It is worthwhile to note that while the smaller kernels achieve higher throughput, they serve as less optimal preconditioners and will require more iterations to converge to the solution of the global problem. In high performance settings, particularly on GPUs, it is 2 WaveHoltz has been demonstrated to be highly effective with implicit time stepping and scaling independently of polynomial degree and mesh scale (see [7]) but we do not consider it here for the complexity it would introduce to the thread-block level code.
typically observed that higher order finite element methods are more performant due to their higher arithmetic intensity. Here however, the DDM cycle has high arithmetic intensity regardless of order, so performance is dictated by asymptotic complexity of the algorithm which is inherently higher for higher order. Observe that the higher order P = 7 discretization attains a peak performance around 200 MDOFs per second. The best performance is achieved by the P = 3 discretization peaking at around 500 MDOFs per second for the range of problems considered. Next, we analyze the cache usage of our method. We measure the memory traffic for the above setup for a 512×512 mesh. Figure 5 shows the share of total memory movement through each level of the memory hierarchy for the B = 1024 kernels. For P = 1 and 3 and T = 1, we see that nearly all of the memory movement takes place in L1 cache with a perfect cache hit rate. For P = 7, T = 1, memory traffic spills over into L2 cache and DRAM. For T = 4 and all P , roughly one third of all data movement is through L1 and half through L2 with the rest spilling out to DRAM. This suggests that the global buffers for T = 4 are too large to sit in cache. This is further emphasized by the cache hit rate with roughly 75% of all L1 fetches missed. The T = 2 kernel finds a balance. Very
P = 1, T = 1
P = 1, T = 4
P = 1, T = 2 0.1%
1.4%
7.8%
32.7%
34.8%
18.5%
65.0% 48.9%
90.8% share of memory traffic cache hit rate L1
98%
L2
L1
74%
L2
75%
98%
P = 3, T = 1
L1
23%
L2
60%
P = 3, T = 2
P = 3, T = 4
0.1%
0.3% 1.8%
32.6%
39.5%
20.0%
60.3% 47.4%
97.9% L1
L1
100%
L2
L1
62%
L2
78%
L2
98%
P = 7, T = 1
26% 63%
P = 7, T = 4
P = 7, T = 2 0.1%
11.1%
14.2%
34.4% 51.6%
L1 L2
37.3%
75% 58%
36.6%
C. Comparison with MINRES
65.5%
L1
49.3%
76%
L2
96%
L1
We compare the runtime of the DDM cycle with five WaveHoltz iterations to the time it takes to evaluate the Helmholtz operator A−iωH −ω 2 M on the entire domain for a constant coefficient problem in [−1, 1]3 with uniform elements and a random input vector. In Figure 6, we plot the relative time of evaluating the DDM cycle compared to evaluating the Helmholtz operator for a representative selection of B and T configurations. The lowest order method (P = 1) is fastest at around 10× operator evaluation time per DDM cycle. The higher order methods are more comparable to two dimensions with all three kernels ranging from 30× to 40× operator evaluation time for the largest problems. In Figure 7, we plot the throughput for increasing problem sizes. The method is roughly half as efficient in three dimensions compared to two. We primarily attribute this to the DOFper-element count scaling as (P + 1)3 in three dimensions, which limits the number of elements per subdomain for a fixed thread block size and reduces the arithmetic intensity of the kernel. Nonetheless, the method remains practical for low-order discretizations in three dimensions.
Cache Level L2
L1
24%
L2
65%
DRAM
Fig. 5. Share of total bytes transferred at each level of the memory hierarchy (L1, L2, DRAM) and cache hit rate for the B = 1024 kernels for p = 1, 3, and 7 in two dimensions.
little memory traffic passes through DRAM, and roughly two thirds of all memory traffic is through L2 with perfect hit rate. These results validate the design principle of the singlethread-block approach. By keeping subdomain data resident in cache, the memory overhead of the computation is substantially reduced. The T = 1 kernel is the ideal configuration when subdomain sizes permit, and T = 2 remains a competitive fallback for larger subdomains. B. Efficiency in Three Dimensions We repeat the two dimensional experiments experiment in three dimensions for meshes ranging from 163 to 1283 elements. In three dimensions, the degrees of freedom per element is (P + 1)3 which implies that our strategy can only be applied to small subdomains. For example, P = 3 requires 64 degrees of freedom per element, so the B = 256, T = 1 kernel configuration can only operate on four elements per subdomain, and the B = 1024, T = 4 configuration only 64 elements. For P = 1, we can operate on subdomains of sizes ranging from 32 elements to 512 elements. Unlike in two dimensions, we do not consider P = 7 which requires 512 threads per element.
Next, we compare the performance of solving the subdomain problem on an A100 with WaveHoltz and MINRES. For a sequence of meshes consisting of n × n elements with P = 3 and we take ω = 0.3n, or roughly ten elements per wavelength. Within each subdomain, both methods iterate until reaching a relative residual of 10−6 in single precision, and 10−12 in double precision. The advantage of WaveHoltz grows with subdomain size for the following reason: as subdomains grow, the local Helmholtz problems become harder and MINRES requires substantially more iterations to converge, with each iteration requiring two inner products and five register arrays that spill into L1 cache. WaveHoltz, by contrast, requires far fewer registers and synchronizations per iteration, and is known to scale well with the problem size [7], so though each iteration is more expensive, far fewer iterations are needed. In addition, MINRES does not benefit meaningfully from mixed precision because the Krylov vectors lose orthogonality, requiring more iterations that nearly cancel out the speedup from faster floating point operations and reduced memory traffic. WaveHoltz, which essentially performs a Fourier transform in time, is numerically stable in reduced precision and benefits directly from the higher throughput. In Figure 8 we plot the runtime for three representative kernels. Across the board WaveHoltz is significantly faster, and the gap widens with subdomain size. Table I quantifies this: in single precision the average speedup ranges from 2× for small subdomains (B = 256, T = 1) to 25× for large ones (B = 1024, T = 2). These speedups are measured in single precision. In double precision the advantage of WaveHoltz is smaller but remains consistent: the speedup ranges from 1.6× to 9× and scales with the subdomain size. The experiment is repeated in three dimensions and the results are summarized in Table I. In three dimensions, the speedup ranges from 0.7× to
B = 256, T = 1
50
DD Time / Operator Time
40
DD Time / Operator Time
30
30
20
30
20
10
20
10
1× 104
105
0
107
106
#DOFs
B = 1024, T = 4
50
40
DD Time / Operator Time
40
0
B = 1024, T = 1
50
10
1× 104
105
#DOFs
Degree P=2
P=1
0
107
106
1× 104
105
106
#DOFs
107
P=4
P=3
Fig. 6. The time to evaluate the DDM cycle (DD time) relative to the time it takes to evaluate A − iωH − ω 2 M on the entire domain (Operator Time) for a sequence of problems of increasing size in three dimensions. P=1
300
250
Throughput [MDOF/s]
Throughput [MDOF/s]
250
200
200
150
150
100
100 50 0
P=3
300
104
105
B = 256, T = 1 B = 256, T = 2
0
106
#DOFs
B = 256, T = 4 B = 512, T = 1
50 106
Kernel B = 512, T = 2 B = 512, T = 4
#DOFs
B = 1024, T = 1 B = 1024, T = 2
107
B = 1024, T = 4
Fig. 7. The throughput (millions of DOFs per second) of the various kernels for fixed degree and increasing problem size. The top panel is p = 1 and the bottom panel is p = 3.
TABLE I AVERAGE SPEEDUP OF THE WAVE H OLTZ SUBDOMAIN SOLVER COMPARED TO THE MINRES SOLVER FOR DIFFERENT KERNEL CONFIGURATIONS IN BOTH SINGLE AND DOUBLE PRECISION IN TWO AND THREE DIMENSIONS . B
104
Runtime [ms]
256
103
512
102
1024 101
105 WaveHoltz, B = 256, T = 1 MINRES, B = 256, T = 1 WaveHoltz, B = 1024, T = 1
106
#DOFs
107
MINRES, B = 1024, T = 1 WaveHoltz, B = 1024, T = 4 MINRES, B = 1024, T = 4
Fig. 8. Runtime comparison of MINRES and WaveHoltz subdomain solvers for various kernel configurations.
T 1 2 4 1 2 4 1 2 4
2D (single) 1.97 7.25 6.52 4.45 11.1 15.4 7.59 24.8 20.1
2D (double) 1.67 2.40 2.56 4.64 3.28 6.60 6.22 8.58 9.10
3D (single) 1.58 2.85 5.03 2.44 3.45 5.87 3.35 4.31 6.48
3D (double) 0.71 2.57 5.91 1.49 3.02 8.48 1.91 4.65 5.18
8.5× with only the (B = 256, T = 1) configuration in double precision resulting in MINRES outperforming WaveHoltz. VIII. C ONCLUSION We presented a procedure for accelerating non-overlapping domain decomposition for the Helmholtz equation on GPUs. Each subdomain is assigned to a single thread block, fusing the entire domain decomposition cycle into a single kernel launch. This replaces expensive global memory traffic with
shared memory and register-level operations, and cache profiling confirms that memory traffic is largely confined to L1 and L2. Second, we use the WaveHoltz iteration as the subdomain solver, which requires a minimal memory footprint and minimal reduction operations, substantially reducing synchronization costs. We considered two threading strategies: one DOF per thread, suitable for subdomains up to 1024 DOFs, and T DOFs per thread (T = 2, 4), which extends the approach to larger subdomains at the cost of increased cache pressure, with T = 2 offering a good balance. Benchmarks on an NVIDIA A100 show that WaveHoltz is 2× to 25× faster than MINRES with the speedup growing as subdomain size increases, and that operating in single precision yields s 2× to 10× speedup over double precision–an advantage largely unavailable to MINRES, whose Krylov vectors quickly lose orthogonality in mixed precision. Several extensions are natural. The zero-order transmission conditions used here can be replaced with optimized Schwarz conditions or higher order ABCs requiring minimal modification to the algorithm and is expected to substantially improve convergence of the outer iteration, particularly in three dimensions where the our method is limited to small subdomains. Further, designing an efficient GPU algorithm for computing and applying coarse spaces by constructing planewave or eigenvector ansatz spaces on the GPU and evaluating the coarse correction within the same kernel framework, for example, is a worthwhile direction and could make the present method robust and competitive for large-scale highfrequency problems. The variable tolerance strategy mentioned in Section V (solving subdomain problems to a weak tolerance initially and tightening as the outer iteration approaches convergence) could substantially reduce total runtime and merits further investigation. Finally, extension to multiple GPUs is natural: the domain decomposition structure already partitions the work into independent subdomain solves, so distributing subdomains across GPUs only requires communication of the interface data with no fundamental change to the single-GPU algorithm.
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
R EFERENCES [1]
O. G. Ernst and M. J. Gander, “Why it is difficult to solve Helmholtz problems with classical iterative methods,” eng, in Numerical Analysis of Multiscale Problems, ser. Lecture Notes in Computational Science and Engineering, Berlin, Heidelberg: Springer Berlin Heidelberg, 2011, pp. 325–363, ISBN: 3642220606. [2] Y. A. Erlangga, “Advances in iterative methods and preconditioners for the Helmholtz equation,” eng, Archives of computational methods in engineering, vol. 15, no. 1, pp. 37–66, 2008, ISSN: 1134-3060. [3] M. Grote and J. Tang, “On controllability methods for the Helmholtz equation,” Journal of Computational and Applied Mathematics, vol. 358, pp. 306–326, 2019.
[12]
[13]
M. J. Grote, F. Nataf, J. H. Tang, and P.-H. Tournier, “Parallel controllability methods for the Helmholtz equation,” Computer Methods in Applied Mechanics and Engineering, vol. 362, p. 112 846, 2020. D. Appelö, F. Garcia, and O. Runborg, “Waveholtz: Iterative solution of the Helmholtz equation via the wave equation,” eng, SIAM journal on scientific computing, vol. 42, no. 4, A1950–A1983, 2020, ISSN: 1064-8275. A. Rotem, O. Runborg, and D. Appelö, “Convergence of the semi-discrete waveholtz iteration,” Journal of Computational Physics, p. 114 882, 2026, ISSN: 00219991. DOI: https://doi.org/10.1016/j.jcp.2026.114882 [Online]. Available: https : / / www. sciencedirect . com / science/article/pii/S0021999126002329 D. Appelö, J. W. Banks, W. D. Henshaw, and D. W. Schwendeman, “An optimal O(N) Helmholtz solver for complex geometry using WaveHoltz and overset grids,” ArXiv, vol. abs/2504.03074, 2025. [Online]. Available: https://api.semanticscholar.org/CorpusID:277596478 B. Deprés, “Méthodes de décomposition de demains pour les problèms de propagation d‘ondes en régime harmonique,” PhD thesis, 1991. M. Taus, L. Zepeda-Núñez, R. J. Hewett, and L. Demanet, “L-sweeps: A scalable, parallel preconditioner for the high-frequency helmholtz equation,” Journal of Computational Physics, vol. 420, p. 109 706, 2020, ISSN: 0021-9991. DOI: https : / / doi . org / 10 . 1016 / j . jcp . 2020 . 109706 [Online]. Available: https : / / www . sciencedirect . com / science / article / pii / S0021999120304800 C. C. Stolk, “A rapidly converging domain decomposition method for the helmholtz equation,” Journal of Computational Physics, vol. 241, pp. 240–252, 2013, ISSN: 0021-9991. DOI: https : / / doi . org / 10 . 1016 / j . jcp . 2013 . 01 . 039 [Online]. Available: https : / / www . sciencedirect . com / science / article / pii / S0021999113000922 A. Royer, C. Geuzaine, E. Béchet, and A. Modave, “A non-overlapping domain decomposition method with perfectly matched layer transmission conditions for the helmholtz equation,” Computer Methods in Applied Mechanics and Engineering, vol. 395, p. 115 006, 2022, ISSN: 0045-7825. DOI: https : / / doi . org / 10 . 1016 / j . cma . 2022 . 115006 [Online]. Available: https : / / www . sciencedirect . com / science / article / pii / S0045782522002420 M. J. Gander and H. Zhang, “A class of iterative solvers for the helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized schwarz methods,” SIAM Review, vol. 61, no. 1, pp. 3–76, 2019. DOI: 10. 1137 / 16M109781X eprint: https : / / doi . org / 10 . 1137 / 16M109781X. [Online]. Available: https://doi.org/10. 1137/16M109781X F. Nataf, “Absorbing boundary conditions and perfectly matched layers in wave propagation problems,”
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
in Direct and Inverse Problems in Wave Propagation and Applications, I. Graham, U. Langer, J. Melenk, and M. Sini, Eds. Berlin, Boston: De Gruyter, 2013, pp. 219–232, ISBN: 9783110282283. DOI: doi:10.1515/ 9783110282283 . 219 [Online]. Available: https : / / doi . org/10.1515/9783110282283.219 T. Hagstrom, R. Tewarson, and A. Jazcilevich, “Numerical experiments on a domain decomposition algorithm for nonlinear elliptic boundary value problems,” Applied Mathematics Letters, vol. 1, no. 3, pp. 299–302, 1988, ISSN: 0893-9659. DOI: https : / / doi . org / 10 . 1016 / 0893 - 9659(88 ) 90097 - 3 [Online]. Available: https : / / www . sciencedirect . com / science / article / pii / 0893965988900973 M. J. Gander, F. Magoulès, and F. Nataf, “Optimized schwarz methods without overlap for the helmholtz equation,” SIAM Journal on Scientific Computing, vol. 24, no. 1, pp. 38–60, 2002. DOI: 10 . 1137 / S1064827501387012 eprint: https://doi.org/10.1137/ S1064827501387012. [Online]. Available: https://doi. org/10.1137/S1064827501387012 P. Chevalier and F. Nataf, “An optimized order 2 (oo2) method for the helmholtz equation,” Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, vol. 326, no. 6, pp. 769–774, 1998, ISSN: 0764-4442. DOI: https://doi.org/10.1016/S0764-4442(98)80047-5 [Online]. Available: https : / / www. sciencedirect . com / science/article/pii/S0764444298800475 Y. Boubendir, X. Antoine, and C. Geuzaine, “A quasioptimal non-overlapping domain decomposition algorithm for the helmholtz equation,” Journal of Computational Physics, vol. 231, no. 2, pp. 262–280, 2012, ISSN: 0021-9991. DOI: https://doi.org/10.1016/j.jcp.2011.08. 007 [Online]. Available: https://www.sciencedirect.com/ science/article/pii/S0021999111004797 R. B. Morgan, “Gmres with deflated restarting,” SIAM Journal on Scientific Computing, vol. 24, no. 1, pp. 20– 37, 2002. DOI: 10.1137/S1064827599364659 M. L. Parks, E. de Sturler, G. Mackey, D. D. Johnson, and S. Maiti, “Recycling krylov subspaces for sequences of linear systems,” SIAM Journal on Scientific Computing, vol. 28, no. 5, pp. 1651–1674, 2006. DOI: 10.1137/040607277 P. Jolivet and P.-H. Tournier, “Block iterative methods and recycling for improved scalability of linear solvers,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’16, Salt Lake City, Utah: IEEE Press, 2016, ISBN: 9781467388153. L. Conen, V. Dolean, R. Krause, and F. Nataf, “A coarse space for heterogeneous helmholtz problems based on the dirichlet-to-neumann operator,” Journal of Computational and Applied Mathematics, vol. 271, pp. 83–99, 2014, ISSN: 0377-0427. DOI: https : / / doi . org/10.1016/j.cam.2014.03.031 [Online]. Available:
[22]
[23]
[24]
[25] [26]
[27]
[28]
[29]
https : / / www . sciencedirect . com / science / article / pii / S0377042714001800 R. Aubry, S. Dey, and R. Löhner, “Iterative solution applied to the helmholtz equation: Complex deflation on unstructured grids,” Computer Methods in Applied Mechanics and Engineering, vol. 241-244, pp. 155– 171, 2012, ISSN: 0045-7825. DOI: https : / / doi . org / 10 . 1016 / j . cma . 2012 . 06 . 007 [Online]. Available: https : / / www . sciencedirect . com / science / article / pii / S0045782512001958 J.-H. Kimn and M. Sarkis, “Restricted overlapping balancing domain decomposition methods and restricted coarse problems for the helmholtz problem,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 8, pp. 1507–1514, 2007, Domain Decomposition Methods: recent advances and new challenges in engineering, ISSN: 0045-7825. DOI: https : / / doi . org/10.1016/j.cma.2006.03.016 [Online]. Available: https : / / www . sciencedirect . com / science / article / pii / S0045782506002696 P. Jolivet, F. Hecht, F. Nataf, and C. Prud’homme, “Scalable domain decomposition preconditioners for heterogeneous elliptic problems,” in Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, ser. SC ’13, Denver, Colorado: Association for Computing Machinery, 2013, ISBN: 9781450323789. DOI: 10 . 1145 / 2503210.2503212 [Online]. Available: https://doi.org/ 10.1145/2503210.2503212 S. G. Johnson, Notes on perfectly matched layers (pmls), 2021. arXiv: 2108.05348 [cs.CE]. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994, ISSN: 0021-9991. DOI: https : / / doi . org / 10 . 1006 / jcph . 1994 . 1159 [Online]. Available: https : / / www . sciencedirect . com / science / article / pii / S0021999184711594 V. Dolean, P. Jolivet, and F. Nataf, “Chapter 2: Optimized schwarz methods,” in An Introduction to Domain Decomposition Methods. 2015, pp. 35–75. DOI: 10 . 1137/1.9781611974065.ch2 eprint: https://epubs.siam. org/doi/pdf/10.1137/1.9781611974065.ch2. [Online]. Available: https://epubs.siam.org/doi/abs/10.1137/1. 9781611974065.ch2 D. A. Kopriva, Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers (Scientific Computation). Springer Netherlands, 2009. K. Ljungkvist, “Matrix-free finite-element computations on graphics processors with adaptively refined unstructured meshes,” in Proceedings of the 25th High Performance Computing Symposium, ser. HPC ’17, Virginia Beach, Virginia: Society for Computer Simulation International, 2017, ISBN: 9781510838222.
[30]
[31]
[32]
[33]
[34]
[35]
[36]
[37]
A. Abdelfattah et al., “Gpu algorithms for efficient exascale discretizations,” Parallel Computing, vol. 108, p. 102 841, 2021, ISSN: 0167-8191. DOI: https : / / doi . org/10.1016/j.parco.2021.102841 [Online]. Available: https : / / www . sciencedirect . com / science / article / pii / S0167819121000879 P. Ghysels and R. Synk, “High performance sparse multifrontal solvers on modern gpus,” Parallel Computing, vol. 110, p. 102 897, 2022, ISSN: 0167-8191. DOI: https: / / doi . org / 10 . 1016 / j . parco . 2022 . 102897 [Online]. Available: https : / / www . sciencedirect . com / science / article/pii/S0167819122000059 NVIDIA, NVIDIA cuDSS (preview): GPU-accelerated direct sparse solver, https://developer.nvidia.com/cudss, Accessed: 2026-03-27, 2025. Y. Saad, “6. krylov subspace methods, part i,” in Iterative Methods for Sparse Linear Systems. 2003, pp. 151–216. DOI: 10 . 1137 / 1 . 9780898718003 . ch6 eprint: https : / / epubs . siam . org / doi / pdf / 10 . 1137 / 1 . 9780898718003.ch6. [Online]. Available: https://epubs. siam.org/doi/abs/10.1137/1.9780898718003.ch6 H. Abdelkhalik, Y. Arafa, N. Santhi, and A.-H. A. Badawy, “Demystifying the nvidia ampere architecture through microbenchmarking and instruction-level analysis,” 2022 IEEE High Performance Extreme Computing Conference (HPEC), pp. 1–8, 2022. [Online]. Available: https : / / api . semanticscholar. org / CorpusID : 251765382 M. Kawai and K. Nakajima, “Low/adaptive precision computation in preconditioned iterative solvers for illconditioned problems,” in International Conference on High Performance Computing in Asia-Pacific Region, ser. HPCAsia ’22, Virtual Event, Japan: Association for Computing Machinery, 2022, pp. 30–40, ISBN: 9781450384988. DOI: 10 . 1145 / 3492805 . 3492813 [Online]. Available: https://doi.org/10.1145/3492805. 3492813 N. Tian, S. Huang, and X. Xu, “Mixed precision blockjacobi preconditioner: Algorithms, performance evaluation and feature analysis,” CCF Transactions on High Performance Computing, vol. 7, no. 2, pp. 114–128, 2025. Y. Guo, E. de Sturler, and T. Warburton, An adaptive mixed precision and dynamically scaled preconditioned conjugate gradient algorithm, 2025. arXiv: 2505.04155 [math.NA]. [Online]. Available: https://arxiv.org/abs/ 2505.04155