arXiv:2606.12850v1 [cs.DC] 11 Jun 2026
High-Order Spectral Element Methods for Wave Propagation on ARM Multicore CPU with SME: Optimizations and Implications Yinuo Wang
Lin Gan
Tsinghua University Beijing, China [email protected]
Tsinghua University Hetao Institute of Mathematics and Interdisciplinary Sciences, Shenzhen Beijing, China [email protected]
Tianqi Mao
Wubing Wan
Zekun Yin
Tsinghua University Beijing, China [email protected]
Tsinghua University Beijing, China [email protected]
Shandong University Jinan, Shandong [email protected]
Wenqiang Wang
Wei Xue
Guangwen Yang
National Supercomputing Center in Shenzhen Shenzhen, China [email protected]
Tsinghua University Beijing, China [email protected]
Tsinghua University Beijing, China [email protected]
Abstract—Wave propagation based on the spectral element method (SEM) is a representative HPC workload, but existing SEM implementations are not well matched to emerging ARM multicore CPUs with Scalable Matrix Extension (SME). We present an SME-enabled optimization of SPECFEM3D on the emerging LX2 processor that combines an SME-aware batched small-matrix kernel for SEM tensor-product operators, a memory-aware hybrid MPI+OpenMP execution scheme for limited-HBM systems, and a dispersion-based iso-accuracy study of the (h, p) tradeoff. At fixed polynomial order, the optimized implementation improves full-application performance by 4–6× over the original code and delivers clear gains over optimized non-SME CPU baselines. Beyond these implementation-level gains, our results suggest that SME shifts the performancefavorable operating point toward higher polynomial orders along the dispersion-based iso-accuracy frontier, further reducing timeto-solution and working-set size. These results indicate that SME affects not only kernel efficiency, but also the practical discretization tradeoff for SEM on modern ARM multicore platforms. Index Terms—Spectral Element Method; Graph Coloring; ARM SME; ARM SVE; Wave Propagation; SPECFEM3D
I. I NTRODUCTION Wave propagation plays a central role in earth modeling, earthquake hazard assessment [1]–[3], and subsurface resource exploration [4], [5]. Because of both its scientific importance and its exceptional demands on computational power and parallelism, numerical wave-propagation simulation has long been a representative large-scale workload in high-performance computing. The pursuit of accurate and efficient modeling of wave propagation in complex media and geometries dates back
to the early development of supercomputers and remains an active challenge today. To meet the need for high accuracy and efficiency in the presence of complex media, irregular geometries, and possible multi-physics coupling, the spectral element method (SEM) has become one of the most promising numerical approaches [6]–[8]. SEM combines high-order polynomial basis functions with tensor-product Gauss-Legendre-Lobatto interpolation on polytope elements, achieving high numerical accuracy together with strong computational efficiency. Compared with classical low-order finite element methods, SEM provides superior accuracy per degree of freedom while retaining the geometric flexibility needed to handle complex domains, boundary conditions, and coupled physics. From an HPC perspective, SEM offers several intrinsic advantages for large-scale simulation. By employing tensorproduct nodal bases with Gauss-Legendre-Lobatto quadrature, it yields a diagonal mass matrix, which eliminates the need for costly mass-matrix inversion and enables efficient explicit time integration. In addition, because only shared nodes on partition boundaries require data exchange, communication is restricted to nearest neighbors, which helps preserve scalability in distributed-memory execution. SEM also benefits from a matrix-free formulation [9], [10], which avoids low-efficiency sparse matrix operations and instead evaluates element contributions on the fly. This reduces memory traffic and better matches the arithmetic structure of high-order methods. As a result, SEM maps naturally to modern SIMD/SIMT architectures and can achieve higher
hardware utilization than traditional sparse finite-element formulations. These properties have made SEM the foundation of several influential high-performance simulation frameworks and applications, including SPECFEM3D [11], MFEM [12], and SEM3D [13], and have contributed to its strong reputation in large-scale scientific computing [14], [15]. Historically, most numerical frameworks were developed under older hardware assumptions or optimized primarily for accelerator-centric environments such as GPU cluster. In contrast, HPC is now undergoing a renewed architectural shift toward large-scale ARM multicore CPU platforms equipped with advanced vector and matrix extensions. For example, China is advancing its national HPC infrastructure with the next-generation ARM-based exascale-class system Lingsheng at the National Supercomputing Center in Shenzhen. Together with established ARM-based flagship systems such as Japan’s Fugaku [16] and its planned successor FugakuNEXT [17], this trend indicates that ARM multicore CPUs will play an increasingly important role in future leading HPC systems. In particular, the ARM Scalable Matrix Extension (SME) introduces matrix-oriented instructions based on outer-product computation and tile storage, providing dedicated hardware support for matrix multiplication [18]. From a hardware perspective, SEM appears to be a natural target for SME, since its stiffness operator is dominated by batched small matrix multiplications arising from tensor-product differentiation and the weak-form operator application. However, this opportunity is not realized directly in conventional SEM configurations. For the polynomial orders commonly used in practice, the resulting matrix sizes are too small to utilize SME tiles efficiently, and therefore a naive SME implementation is often unable to deliver a substantial advantage over conventional SIMD-based kernels. To address this challenge, we develop an SME-enabled batched small matrix multiplication kernel based on software pipelining and vector aggregation, which improves SME tile utilization and achieves higher performance than traditional SIMD-backed implementations. More importantly, once this kernel-level inefficiency is removed, the limiting factor is no longer only the microarchitectural mapping, but also the application-level discretization choice. For wave-propagation problems, different combinations of element size h and polynomial order p form an iso-accuracy frontier w.r.t. dispersion error. Since SME changes the efficiency of the dominant batched matrix multiplication kernels along this frontier, we observe that its introduction shifts the performancefavorable operating point toward the higher-p side. We validate this observation through dispersion-based iso-accuracy wave-propagation experiments in homogeneous media.In the tested equal-accuracy regime, higher-p discretizations paired with SME-enabled kernels reduce both time-to-solution and working-set size. To support realistic wave-propagation simulations, another major challenge must be addressed: memory capacity. Beyond the distributed field variables, practical SEM applications such
as SPECFEM3D typically require each rank to maintain substantial auxiliary data, including tomographic medium descriptions, surface geometry information, preprocessing metadata, and communication buffers. Under the traditional MPI-only execution model adopted by many existing SEM frameworks, these per-rank data structures are replicated across a large number of processes. On modern ARM multicore platforms with very high core counts, such replication places severe pressure even on the DDR memory subsystem. During the PDE solution stage, additional memory is consumed by halo buffers, duplicated boundary degrees of freedom, and processlocal runtime metadata, further reducing the limited HBM capacity available on each NUMA domain. We therefore develop a hybrid MPI+OpenMP parallelization scheme for SEM on multicore CPUs. MPI processes are mapped to NUMA domains, while thread-level parallelism is exploited within each domain to reduce per-process memory replication. To eliminate shared-memory race conditions during nodal updates, we apply graph coloring instead of relying on costly atomic operations. At the inter-process level, we introduce a dedicated communication thread to achieve fully asynchronous neighbor exchanges and overlap communication with computation. This hybrid execution scheme substantially reduces memory consumption and enables realistic wavepropagation simulations on high-performance multicore CPUs. In this work, we select SPECFEM3D as a representative wave-propagation application and implement our proposed optimizations on LX2, the processor powering China’s new supercomputer Lingsheng. Our optimized implementation achieves a full-application performance improvement of 4–6× over the baseline, while enabling realistic wave-propagation workloads under the memory constraints of the target ARM multicore platform. Our main contributions are summarized as follows: 1) We design an SME-aware batched small-matrix kernel for SEM tensor-product operators on ARM multicore CPUs. By combining software pipelining with layoutaware vector aggregation, the proposed kernel improves utilization of SME tiles for the practically relevant small operator sizes arising in SEM wave propagation problem. 2) We show that SME changes the practical efficiency trend of SEM across polynomial orders, making higherp discretizations substantially more attractive on CPUs, and experimentally use a dispersion-based iso-accuracy study to show that SME shifts the hardware cost structure of the dominant tensor-product operators and thus motivates renewed consideration of higher-order SEM on modern ARM CPUs. 3) We develop a memory-aware hybrid MPI+OpenMP execution scheme for SPECFEM3D on multicore NUMA systems. This design reduces per-process memory replication and supports realistic wave-propagation workloads under limited HBM capacity, using graph coloring to remove shared-memory update conflicts and a dedicated communication thread to overlap neighbor
exchange with computation. η-axis
ζ-axis
II. BACKGROUND A. LX2 Architecture
ξ-axis
Die 0 MCM IO
IO HBM
Die 1
UC
...
UC
...
IO
IO
HBM
HBM
UC
HBM
UC
UC
...
UC
...
UC
...
HBM
HBM
UC
... ...
NUMA 1
NUMA 5
...
...
UC
...
HBM
HBM
UC
UC
HBM
UC
HBM
...
UC
HBM
UC
HBM
NUMA 6
NUMA 2
HBM
...
IO
NUMA 4
NUMA 0
HBM
Evaluation of the ξ-Direction Derivative
MCM IO
...
UC
...
HBM
HBM
UC
... ...
NUMA 3
NUMA 7
DDRC
DDRC
Fig. 1. LX2 CPU Architecture
The LX2 is a dual-die ARMv9-A processor featuring 304 physical cores operating at 1.6 GHz. Each die comprises four NUMA domains of 38 cores. Every core includes private L1/L2 caches and supports the Scalable Matrix Extension (SME) with specialized outer-product units. The memory subsystem features 4 GB of High Bandwidth Memory (HBM) per NUMA node (500 GB/s), configurable as either a hardware-managed cache or a software-addressable flat memory. Each die also accesses 1 TB of shared DDR memory (120 GB/s). A 160-channel System Direct Memory Access (SDMA) engine facilitates asynchronous transfers between DDR and HBM, effectively overlapping memory traffic with computation. For intensive matrix computations, the SME unit employs a vector outer-product model. It computes the outer product of two vectors loaded from SVE registers, accumulating the result into a 64 × 64-byte 2D architectural register known as the ZA tile. The unit maintains full IEEE floating-point compliance. Data moves flexibly between the ZA tile and SVE registers via horizontal or vertical slice instructions. By interleaving operations to hide instruction latency, SME delivers approximately 4× the theoretical peak performance of standard SVE Multiply-Accumulate (MLA) instructions. B. Spectral Element Method The spectral element method (SEM) is a class of highorder finite element methods. On a three-dimensional mesh it typically employs hexahedral elements. Within each element on the reference domain with coordinates (ξ, η, ζ), the interpolation (nodal) points are defined as the tensor product of onedimensional Gauss-Legendre-Lobatto (GLL) points. If {ξi }N i=0 denotes the set of (N + 1) one-dimensional GLL nodes, then the coordinates of the three-dimensional interpolation nodes in reference space are xijk = ξi , ξj , ξk , for i, j, k = 0, . . . , N .
Fig. 2. ξ-Derivative in Tensor-Product Element Basis
The corresponding nodal basis functions within each element are given by the tensor product of one-dimensional Lagrange interpolation polynomials: lijk (ξ, η, ζ) = li (ξ) lj (η) lk (ζ), where li (ξ) is the i-th one-dimensional Lagrange basis polynomial associated with the GLL nodes {ξ0 , ξ1 , . . . , ξN } for a polynomial of degree N . It can be written explicitly as li (ξ) =
(ξ − ξ0 ) · · · (ξ − ξi−1 )(ξ − ξi+1 ) · · · (ξ − ξN ) . (ξi − ξ0 ) · · · (ξi − ξi−1 )(ξi − ξi+1 ) · · · (ξi − ξN )
By construction, these basis functions satisfy the Kronecker delta property li (ξj ) = δij ,
i, j = 0, . . . , N,
which implies lijk (ξp , ηq , ζr ) = δip δjq δkr for the threedimensional tensor-product basis. This tensor-product interpolation basis has a particularly convenient structure. For a function with nodal values fijk = f (ξi , ξj , ξk ), its value at an arbitrary point (ξ, η, ζ) in the reference element can be written as f (ξ, η, ζ) =
N X
fijk li (ξ) lj (η) lk (ζ).
i,j,k=0
Taking derivatives in ξ-direction with respect to the reference coordinates yields N X ∂f (ξ, η, ζ) = fijk li′ (ξ) lj (η) lk (ζ) ∂ξ i,j,k=0
At the interpolation points (ξp , ξq , ξr ), the Kronecker-delta property lj (ξq ) = δjq and lk (ξr ) = δkr implies that, in each derivative, only nodes aligned with ξ-axis contribute to the result, while nodes with mismatched indices have no contribution. Consequently, evaluation of the ξ-direction derivative can be expressed as a sequence of small matrixmatrix multiplications of an (N +1)×(N +1) one-dimensional differentiation matrix with an (N +1)×(N +1) slice of nodal values taken along the ξ-direction, and similarly for η and ζ (see Fig. 2).
For the representative application SPECFEM3D we consider, the governing equation is the elastic wave equation
for all test functions w. In the SEM discretization, both the trial and test functions are chosen as the tensor-product Lagrange basis functions described above. The left-hand side gives rise to the mass matrix; with inexact numerical integration at the Gauss-LobattoLegendre nodes, this matrix is diagonal to machine precision. The main computational hot spot is the first term on the righthand side, which corresponds to the stiffness matrix. SEM employs a matrix-free approach to evaluate the action of this stiffness operator. The computation of the stress tensor T involves spatial derivatives of s in the three reference directions, which, due to the tensor-product structure, can be expressed as three small batched matrix-matrix multiplications. Applying the weak form, i.e., contracting with the gradients of the test functions, leads to an additional three small batched matrixmatrix multiplications of the same structure.
Gordon Bell Prize-winning real-time tsunami forecasting application, which employed polynomial spaces up to fifth order [15]. In contrast, high-order tensor-product optimization on CPUs has received comparatively less attention, with many implementations continuing to rely on general-purpose densekernel libraries, such as LIBXSMM [21]. Graph coloring has long been used to eliminate write conflicts during shared-memory assembly and element updates, and communication overlap via thread-based progress has been extensively studied in MPI runtimes and applications [22]– [25]. While both graph coloring and progress threads are not novel concepts, our work integrates them into a unified SEM execution scheme for realistic wave-propagation workloads on ARM multicore CPUs, where CPU-side atomic updates are costly. Additionally, we employ a dedicated communication thread to ensure asynchronous progress under a hybrid MPI+OpenMP model [26]–[28]. In comparison to prior high-order FEM/SEM optimization work, our contribution is twofold. First, we focus on ARM multicore CPUs with SME, a setting that has received less attention than GPU-oriented SEM optimizations. Second, rather than solely optimizing the tensor-product kernel, we demonstrate that SME shifts the performance-favorable operating point along the dispersion-based iso-accuracy (h, p) frontier for wave-propagation SEM. Thus, our work is not merely a microarchitectural kernel optimization but an application- and architecture-aware co-design of discretization choices, kernel mapping, and hybrid parallel execution.
III. R ELATED W ORK
IV. O PTIMIZATION
A diverse range of frameworks exists for spectral element methods (SEM) in wave propagation, including specialized packages such as SEM3D [13] and SPECFEM3D [11], as well as more general high-order finite-element libraries such as MFEM [12] and DEAL .II [19], which prioritize generality and performance portability. Among these, SPECFEM3D stands out as one of the most influential frameworks for wave propagation. This widely used, open-source seismic simulation tool, based on the spectral-element method, is designed for large-scale applications, ranging from local and regional models to global Earth simulations. Its significance in both the seismology and HPC communities is underscored by its receipt of the ACM Gordon Bell Award in 2003, its selection as a Gordon Bell finalist in 2008, and its distribution through the Computational Infrastructure for Geodynamics (CIG) [14]. More broadly, high-order tensor-product FEM/SEM has garnered sustained attention in the HPC community due to its dense operator structure and matrix-free formulation, which align well with modern hardware architectures. Recent efforts, such as LIB CEED [20] and MFEM, have focused on matrixfree and partially assembled high-order discretizations as core optimization strategies, particularly for GPU platforms, with LIB CEED serving as a tensor-product backend for MFEM’s high-order operator implementation. These efforts have shown notable application-level impact: MFEM formed part of the high-order finite-element infrastructure for the 2025 ACM
A. Aggregated SME Kernels for Batched Small Matrix Multiplications
ρ ∂t2 s = ∇ · T + f ,
T = C : ∇s,
where ρ is the density of the solid, s is the displacement field, T is the stress tensor governed by Hooke’s law, and C is the fourth-order stiffness tensor. The corresponding weak form reads Z Z Z ρ w · ∂t2 s d3 x = − ∇w : T d3 x + w · f d3 x, Ω
Ω
Ω
In this section, we address the mismatch between the SME ZA tile size and the dimensions of the batched small matrix multiplications in SEM kernels by combining software pipelining with SVE aggregation scheme. After the global degrees of freedom are gathered, each field component is stored in an element-local (N +1) × (N +1) × (N +1) cube with the ξ-direction as the leading dimension (optionally with padding), and the one-dimensional differentiation operator is stored as an (N +1) × (N +1) coefficient matrix. We target polynomial orders from p = 4, a common setting in existing SEM codes, up to p = 15, which fully matches the ZA tile size for single-precision operands. Among the three reference directions, the η-direction derivative most closely resembles a batched matrix-matrix multiplication, where each ζ-slice of the element-local cube multiplies the one-dimensional coefficient matrix. SME exposes multiple ZA tiles that can be used for software pipelining to hide the latency of outer-product instructions. For large matrix multiplications, the four ZA tiles are typically grouped as a 2 × 2 tile set for a single matrix pair, but this is unnecessary for batched small matrix multiplications whose dimensions are at most the vector length. Instead, we map different ζslices to different ZA tiles and interleave outer products across the batch. Each tile processes a distinct slice, so all tiles
η-Direction Derivative
ξ-Direction Derivative
Data Layout in Memory
Nodal Values Slices
For N+1 ≤ VL/2, one SIMD register can hold two η-lines
Load Vector η-th Line
Load Vector
(η+1)-th Line
(ζ, η)
(ζ, η+1)
(ζ+1, η)
(ζ+1, η+1)
(ζ+1)-th Slice
00
01
02
03
00
01
02
03
10
11
12
13
10
11
12
13
20
21
22
23
20
21
22
23
30
31
32
33
30
31
32
33
In Memory
(ζ, η)
(ζ+1, η)
(ζ, η+1)
(ζ+1, η+1)
Transpose Required Vector Discontinuous
ζ-Direction Derivative Load Vector Extract Slice along η-axis
Derivatives
Horizontal Load
TLB shuffle
ζ-th Slice
Differentiation Matrix
ZA Tile 00
10
20
30
01
11
21
32
02
12
22
32
03
13
23
33
In Memory
Naturallly Concatenated (ζ, η)
(ζ, η+1)
Data Layout in Memory
Fig. 3. Aggregated SME-enabled Batched Matrix Multiplication
remain active and outer-product latency is hidden by switching between slices, improving instruction throughput for the SEM workload. We further increase ZA utilization within each ζ-slice for low polynomial orders. Algebraically, this is equivalent to concatenating multiple ζ-slices along the ξ-axis and applying the same coefficient matrix to the concatenated data. We realize this logical concatenation on the fly using the vectorshuffling capabilities of the Scalable Vector Extension (SVE). For N = 4-7, a single SIMD vector register can hold two ηlines. A vector load starting at one η-line brings both that line and the subsequent η+1 line. By combining this vector with the corresponding vector from the next ζ-slice and applying an SVE table-lookup shuffle, we construct operands for two consecutive outer-product iterations, so a single outer-product updates two ζ-slices and effectively doubles ZA utilization. The ξ-direction derivative has the same algebraic structure as the η-direction case: each ζ-slice is multiplied by the same coefficient matrix. To improve SME utilization, we aggregate all ζ-slices along the first dimension, forming a larger ((N + 1)2 ) × (N + 1) matrix in row-major order. However, this inmemory layout does not match the operand requirements of SME outer-product instructions, since each column is stored with non-unit stride. To avoid strided loads, we apply an inregister matrix transpose using the SME vertical and horizontal move instructions. This transformation reorganizes the data so that each column becomes contiguous in memory, enabling efficient SME outer-product evaluation. For the ζ-direction derivative, the operation can be viewed as a batched matrix multiplication between the coefficient matrix and a sequence of η-slices. We aggregate these η-slices
along the second dimension, forming an (N + 1) × ((N + 1)2 ) matrix. In this arrangement, each row corresponds to a contiguous ξη-plane in memory, so the resulting matrix multiplication matches the SME operand layout directly and can be executed without additional transpose or shuffle overhead. We next derive a simple attainable-performance upper bound for the SME-enabled kernels based on instructionthroughput considerations. For the ζ- and η-direction derivatives, evaluating all derivatives requires N +1 outer-product operations. Since the matrix dimensions are smaller than the SME tile size, each outer product utilizes only a fraction NV+1 of the full tile capacity. In L addition, N + 1 outer-product updates are followed by N + 1 store operations to write back the results. Under this simplified CPI-based view, the attainable floating-point throughput is upper-bounded by CP IOP A N +1 · F LOP SSM E . CP IOP A + CP ILSU VL For the ξ-direction derivative, the in-memory data layout does not match the operand layout required by SME outerproduct instructions, so an additional transpose is needed. Using a rough throughput estimate, transposing an (N + 1) × (N + 1) matrix introduces another N + 1 loads and N + 1 stores through the ZA tiles. The corresponding attainableperformance upper bound is therefore reduced to CP IOP A N +1 · F LOP SSM E . CP IOP A + 3 CP ILSU VL Although simplified, these expressions make clear that before tile saturation, the maximum attainable SME throughput
scales approximately linearly with NV+1 . This observation moL tivates the next subsection, where we examine how improved efficiency at higher polynomial orders can shift the favorable operating point on the (h, p) design frontier.
this regime, the use of higher polynomial order together with SME-enabled kernels becomes particularly attractive on ARM multicore platforms.
B. SME and the Performance-Favorable Operating Point on the (h, p) Iso-Accuracy Frontier
C. Hybrid SEM Parallelization Scheme for Multicore CPUs
From Sec.IV-A, our SME-enabled kernel maps the SEM batched small matrix multiplications to SME tiles with a and approaches the corresponding utilization factor of NV+1 L theoretical upper bound. Beyond this point, further improvement in SME efficiency cannot be achieved solely at the kernel level, and instead requires revisiting the applicationlevel discretization choice. For wave-propagation problems, the two dominant numerical error sources are temporal dispersion and spatial dispersion. The former arises from discrete time integration, while the latter is caused by insufficient spatial resolution to accurately represent high-frequency wave components. To reduce spatial dispersion, two standard refinement strategies are available: h-refinement, which decreases the element size, and p-refinement, which increases the polynomial order. Accordingly, the set of (h, p) pairs that achieve the same spatialdispersion accuracy defines an iso-accuracy frontier with respect to dispersion error. For smooth solutions, SEM exhibits rapid convergence with increasing polynomial order, so higher p can achieve the same spatial dispersion accuracy with fewer elements. Historically, however, the O(N 4 ) cost of the tensor-product stiffness operator has made high-p configurations less attractive in practice. SME changes this trade-off. Since SME efficiency increases approximately linearly with polynomial order before saturation at the vector length, the effective hardware cost scaling of the dominant kernel is reduced from 4 batched matrix multiplication 3 N N O h3 to approximately O h3 in this regime. As shown on Fig. 4, the performance-favorable operating point along the iso-accuracy frontier shifts toward the higher-p side. Moving toward higher p also brings an additional memory benefit. Because higher-order SEM can achieve the same spatial accuracy with fewer points per wavelength, the total number of global degrees of freedom can be reduced, leading to a smaller working set in memory. For ARM platforms such as LX2, where fast HBM capacity is limited, this directly translates into support for larger simulation domains or fewer required processors. There is, however, an important counter-effect. The minimum nodal spacing, which governs the maximum allowable time step through the CFL condition, decreases as O(N −2 ) for Gauss–Lobatto–Legendre nodes. In principle, this could require more time steps for a fixed physical simulation time and offset the gains from improved SME utilization. As we show in the experimental section using homogeneous equalaccuracy wave-propagation experiments however, for applications with high accuracy requirements such as reverse-time migration, the time-step restriction needed to control temporal dispersion is already more stringent than the CFL bound. In
On modern multicore CPUs, growing core counts make MPI-only execution increasingly inefficient for SEM. In application frameworks such as SPECFEM3D, each process maintains substantial private data, including tomographic models, surface topography, and preprocessing metadata, so launching many MPI ranks creates significant pressure even on DDR memory. This pressure is more severe at the NUMA and HBM levels. On LX2, each NUMA domain provides only 4 GB of HBM. With one MPI process per core, replicated boundary DoFs, communication buffers, and per-process metadata can consume a substantial fraction of this capacity, leaving little space for the PDE solve and, under huge-page mode, even risking allocation failure. These considerations make hybrid parallelization essential for multicore CPUs. In our design, we map MPI processes to NUMA domains and employ multi-threaded parallelism within each NUMA. The whole workflow is shown in Fig. 5. On the NUMA-local, multi-threaded level, race conditions become a central issue, since adjacent elements share global nodes that must be updated consistently. In SPECFEM3D, OpenMP pragmas with atomic updates are used to remove these races. However, unlike GPGPUs, which often provide hardware floating-point atomics, CPUs implement atomics via load-linked/store-conditional (LL/SC) or compare-and-swap loops that repeatedly retry the update if another thread has modified the target location. Such atomic operations disrupt instruction pipelining and expose the full HBM access latency to the application, leading to poor scalability on multicore CPUs. To avoid this overhead, we adopt a graph-coloring approach to eliminate race conditions without atomics. SPECFEM3D uses a conformal hexahedral mesh, so elements can only be connected through shared vertices, edges, or faces. We exploit the existing global-to-local node mapping to construct an element-to-corner-vertex map, and then build an element adjacency graph via a sparse matrix-matrix multiplication on this incidence structure. A graph-coloring algorithm is applied to this adjacency graph, and elements are executed color by color. Within each color class, no two elements share nodes, so threads can update all nodal values without atomics or synchronization, removing shared-memory race conditions while preserving high pipeline and cache efficiency. At the multi-process level, truly asynchronous neighbor exchanges become essential. However, the MPI progress model is often a limiting factor: nonblocking communication is only guaranteed to make progress when the application periodically invokes MPI routines such as MPI_test or MPI_wait. In practice, these calls may block the invoking thread or make only limited internal progress.
(a) Computation Cost on Iso-Accuracy Frontier w/o SME (b) Computation Cost on Iso-Accuracy Frontier with SME 160
160
Inaccurate Region
120
cost=0.0005
120
cost=0.01
cost=0.001
Inaccurate Region
h 80
cost=0.0013
h 80 Feasible Region cost=0.1
40
Feasible Region cost=0.01 40
cost=1 cost=10
04
6
8
p
10
12
04
14
Iso-accuracy frontier
Feasible region
Inaccurate region
cost=0.1 cost=1 6
8
p
10
12
14
Iso-computation-cost line
Fig. cost scales as 4. Shift of the Performance-Favorable Operating Point on the (h, p) Iso-accuracy Frontier. Without SME, the dominant computational p4 p3 O h 3 , as shown in (a). With SME, improved utilization of tensor-product operators reduces the effective cost scaling to O h3 , shifting the favorable operating point toward higher polynomial orders, as shown in (b).
Graph Coloring
Decompose
Domain Partition Communicate Thread
Subdomain Mesh
IDLE
Init Comm
Atomic Notify
Compute Threads
...
Time Line
Phase 1: Outer Elements
...
...
Graph Coloring
...
...
...
Outer Elements Inner Elements (Phase 1) (Phase 2)
Drive MPI Progress
...
Phase 2: Inner Elements
Atomic Notify
...
...
Wait
Next Step...
Fig. 5. Hybrid Multiprocess/Multithread SEM Parallelization Scheme
Given the abundant parallel resources on modern multicore CPUs, we dedicate one CPU core per MPI process to a progress thread responsible for handling all communication. The main thread initiates nonblocking neighbor exchanges and communicates with the progress thread via atomic variables. The progress thread repeatedly drives MPI progress on behalf of the process, while the main thread continues with computation. This design provides genuine asynchronous communication and enables full overlap between communication and computation for the SEM neighbor exchanges. D. Code-Specific Optimizations for SPECFEM3D To further exploit the microarchitecture of the LX2, we refactor the core computing routines of SPECFEM3D. First, for the kernel computing spatial derivatives, strains, and stresses, we fuse its multiple isolated spatial loops into a single grid-point traversal. This structural change eliminates intermediate gradient arrays, using scalars instead to avoid memory bandwidth bottlenecks. Furthermore, we move the element-regularity conditional completely outside the computational loops and specialize the kernel into two separate execution paths. By eliminating control-flow divergence within
the hot loops, this transformation removes an obstacle to compiler vectorization. Second, for the gather and scatter routines, we collapse the deeply nested 3D loops over the Gauss-Lobatto-Legendre points into a single 1D loop. This structural flattening enables higher SIMD utilization. V. E XPERIMENTAL E VALUATION A. Experimental Setup We use the LX2 platform described in Sec.II-A as our primary experimental system. Unless otherwise stated, all experiments are conducted on a single server equipped with two LX2 CPUs. We compile the code with the vendorprovided BiSheng 5.0.0 compiler and use Hybrid MPI 25.2.0. As the application baseline, we use the main branch of SPECFEM3D from GitHub at commit 7a1c764. To ensure fair comparison across optimization stages and baselines, we fix the process/thread mapping, NUMA placement, and memory allocation policy throughout all experiments. For cross-platform comparison, we additionally use a server equipped with two Intel Xeon Gold 6248R processors (48 cores in total), 1.0 TiB DDR4 memory at 3200 MT/s, a 1.0 TiB NVMe SSD, and four NVIDIA A100-PCIE-40GB GPUs. We compile each implementation on different platforms with the highest stable optimization level together with architecturespecific code generation and fast-math enabled. B. Kernel-Level Evaluation of SME-Enabled Batched Small Matrix Multiplication In this subsection, we evaluate the proposed SME-enabled batched small matrix multiplication kernels. We use the effective single-core FLOP rate as the primary performance metric. As baselines, we include both the compiler-optimized implementation and a hand-written highly optimized SVE-intrinsic
(a) -summation
(b) -summation
(c) -summation
Performance (GFLOPS)
250 200 150 100 50 0
4 5 6 7
9
11
Polynomial order N
13
15
4 5 6 7
9
11
Polynomial order N
Compiler-Opt
SVE
SME
13
15
4 5 6 7
9
11
Polynomial order N
13
15
SME Theoretical Bound
Error
Fig. 6. Batched Small Matrix Multiplication Performance
0.0200 0.0175 0.0150 0.0125 0.0100 0.0075 0.0050 0.0025 0.0000
(a) Time step (P=7)
(b) Element size (P=4)
t=5e-5 t=1e-4 t=2e-4 t=3e-4 t=4e-4 t=8e-4
0
5000
(c) Element size (P=7)
h=17.5 h=20 h=22.5 h=25
10000
Propagation distance (m)
0
5000
(d) Element size (P=15)
h=40 h=45 h=50 h=60
10000
Propagation distance (m)
0
5000
h=110 h=120 h=125 h=130 h=135
10000
Propagation distance (m)
0
5000
10000
Propagation distance (m)
Fig. 7. Wave Propagation Error on Different (h, p) Settings
implementation1 . To assess how close the SME kernels are to the architectural limit, we also report the theoretical upper bound derived in Sec.IV-A. As shown in Fig. 6, the ζ-direction kernel does not require transpose or explicit aggregation and therefore achieves the highest performance. For p = 7 and p = 15, where the matrix dimensions align well with the hardware, the achieved performance is close to the theoretical upper bound. In contrast, for the commonly used p = 4 configuration, SME provides little advantage over SVE, which is consistent with the analysis in Sec.IV-A. We also observe that, although the SVE kernel performance increases with p, it still reaches less than 60% of the theoretical SVE compute capability. This is expected, since high SVE throughput requires sustained dual-FMA issue, whereas in small matrix multiplication the setup and store overheads remain significant. For the η-direction kernel, the proposed aggregation scheme allows performance comparable to that of the ζ-direction kernel. In these two directions, the advantage of SME over 1 We additionally evaluated L IB XSMM backend as a reference. However
this generic backend was consistently below our specialized non-SME vectorized implementation due to lack of static shape information, so we report the latter as the strongest CPU baseline for assessing SME-specific gains.
SIMD-based implementations is substantial: for p = 7 and p = 15, the SME kernels achieve up to 4× speedup over the SIMD implementation. The ξ-direction kernel incurs additional transpose overhead and therefore performs worse than the other two directions. Nevertheless, it still outperforms the SIMD-based implementations when p > 6, and achieves approximately 2× speedup at p = 15. C. (h, p) Dispersion-Based Iso-Accuracy Frontier Study 1) Equal-Accuracy Configuration Selection: To quantify dispersion error and construct dispersion-based iso-accuracy configurations, we consider acoustic wave propagation in a homogeneous medium and record the waveform at stations placed along one spatial direction. The simulation error is measured with respect to the analytical solution using the energy norm R 2 (w(x, t) − wtheory (x, t)) dt R . error(x) = wtheory (x, t)2 dt The medium has density 1000 kg/m3 and P-wave speed 1500 m/s. We use a Ricker source with dominant frequency 20 Hz, the acoustic solver in SPECFEM3D, and second-order Newmark time stepping.
Time to solution (s)
Memory (GB)
P = 4 vs P = 7 60
P = 7 vs P = 15 60
56.0 GB
40
40 17.6 GB
20 0
600
P=4
P=7
0 300
536.6 s
200
200
100
139.0 s
P=4
P=7
27.8 GB
20
400 0
52.8 GB
0
P=7
P=15
Elapsed time (s)
(a) P = 4 1500 1000 500 0
1000 800 600 400 200 0
178.1 s
P=15
Original
Vector
Original
Vector
Original
Vector
Hybrid
SME
Hybrid
SME
Hybrid
SME
(b) P = 7
(c) P = 15
1500 1000 500 0
250.8 s
P=7
D. End-to-End Performance of Optimized SPECFEM3D
Elapsed time (s)
In our high-accuracy setting, the temporal-dispersion constraint is more restrictive than the CFL limit. As shown in Fig. 7, the maximum admissible element size is 20 m for p = 4, 50 m for p = 7, and 130 m for p = 15. 2) Computational and Memory Cost on the Iso-Accuracy Frontier: Using the equal-accuracy (h, p) configurations derived above, we evaluate isotropic elastic wave propagation and compare two pairs of discretizations: p = 4 versus p = 7, and p = 7 versus p = 15, with our optimized SPECFEM3D. The goal is to quantify the additional application-level benefit that becomes available once SME-aware execution makes higher-order discretizations practically attractive. In each case, the computational domain is chosen to saturate the available HBM resources under the lower-order configuration, and the element sizes for the higher-order configurations are scaled according to the equal-accuracy ratios determined above. Because the temporal-dispersion constraint is stricter than the CFL stability limit in this regime, we use the same time step for all configurations and run the SPECFEM3D solver for 15,000 steps, recording the total simulation time.
order setting, workloads that cannot fit within a single server node at low polynomial order can be executed within one node. The trend continues when moving from p = 7 to p = 15, yielding a further 1.8× reduction in memory usage and a 1.4× reduction in time-to-solution. These results support the use of higher polynomial order together with SME-enabled kernels for wave-propagation simulations on ARM multicore platforms.
Elapsed time (s)
Motivated by the accuracy requirements of reverse-time migration, we require that after 8 s of propagation the waveform error remain below 1%. Since wave-propagation simulations are affected by both temporal and spatial dispersion, we first determine the maximum allowable time step by running on a sufficiently fine mesh such that spatial dispersion is negligible. As shown in Fig. 7, satisfying the target accuracy requires ∆t = 1 × 10−4 s. We then determine the admissible element size for p = 4, 7, 15 by testing different element sizes and using ∆t = min ∆tCFL , 1 × 10−4 .
Original Stiffness
Communication
SME (w/o color) Mass + Update
Fig. 8. SME-enabled Performance Comparison on Iso-Accuracy Frontier
Fig. 9. Performance Breakdown
As shown in Fig. 8, switching from p = 4 to p = 7 yields a 3.18× reduction in memory usage and a 3.85× reduction in time-to-solution. The memory reduction arises because higher polynomial order can achieve the same accuracy with fewer points per wavelength, thereby reducing the total number of global degrees of freedom. This is particularly important for applications such as reverse-time migration, where wavefield history must be stored, and for platforms such as ARM multicore CPUs with limited HBM capacity. Under this higher-
1) Breakdown of Optimization Contributions: We evaluate the contribution of each optimization to reducing time-tosolution under three polynomial-order settings: p = 4, 7, 15. The earlier stages quantify the benefit of general CPU-side optimization and hybrid execution at fixed discretization, whereas the final stage isolates the incremental contribution of SME relative to previous stages. For each setting, we choose a mesh that saturates the available HBM capacity on one server node and run the simulation for 15,000 steps.
Performance (FLOPS)
8
Performance (FLOPS)
(a) Single-node Workset Scaling
1e12
7 6 5 4 3 2
10
4.0 1e13 3.5 3.0 2.5 2.0 1.5 1.0 0.5 1
Performance (FLOPS)
The vectorization optimization in Sec.IV-D improves performance by 2.1×, 2.3×, and 1.4× over the original SPECFEM3D for p = 4, p = 7, and p = 15, respectively, and thereby establishes a stronger non-SME baseline for the subsequent stages. Graph coloring further improves performance by 1.40× and 1.47× for p = 4 and p = 7, but yields only 0.94× for p = 15. The degradation at p = 15 can be explained by two factors. First, when the batched small matrix multiplication becomes sufficiently expensive, the cost of serialized atomic updates accounts for a smaller fraction of the total runtime. Second, executing the stiffness operator color by color reduces the effective working set of each OpenMP region and weakens the locality benefits introduced by DoF renumbering. The final stage adds SME on top of the optimized nonSME implementation. For the stiffness operator alone, SME yields speedups of 1.06×, 1.27×, and 7.69× for p = 4, p = 7, and p = 15, respectively; these translate into endto-end simulation speedups of 1.05×, 1.26×, and 5.92×. The improvements at p = 4 and p = 7 are consistent with the kernel-level results. For p = 15, the observed application-level speedup is larger than the kernel-level SME gain alone would suggest. The reason is that the baseline implementation uses a generic runtime-parameterized batched matrix multiplication path for polynomial orders higher than 10, whereas our optimized implementation combines SME-specific acceleration with static shape specialization. Accordingly, the improvement at p = 15 reflects the joint benefit of architectural support and specialization to SEM operator structure. The communication thread provides a mechanism to overlap communication with computation. In practice, however, its impact is limited. Domain decompositions produced by METIS or Scotch typically introduce a workload imbalance of about 5–10%, and under the hybrid parallelization scheme each process exchanges less than 1 MB of halo data compared to a working set of several gigabytes. Consequently, the remaining asynchronous wait time is dominated by load imbalance rather than communication overhead, and the communication-thread optimization contributes little to end-to-end performance. Finally, to confirm the effectiveness of graph coloring at p = 15, we conduct an ablation study on the final optimized version by removing graph coloring and reverting to OpenMP atomic updates. This increases the simulation time by 1.16×, demonstrating that graph coloring remains beneficial once the kernel is optimized. 2) Scaling Experiments: For the intra-node scaling experiments, we additionally include a cross-platform comparison against two NVIDIA A100 GPUs using the same mesh configuration. To present results for different polynomial orders in a single figure, we use the memory footprint on a LX2 server node as the horizontal axis. As shown in Fig. 10, higher polynomial order increasingly favors the ARM platform because SME utilization improves with p. For p = 4, the A100 is approximately twice as fast, reflecting the low SME utilization in this regime. At p = 7, LX2 becomes competitive with the A100. At p = 15, the A100 implementation can no longer be compiled because
1.6 1e14 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0 1
20
30
40
50
Memory consumption (GB)
(b) Strong Scaling 40% 61%
78% 91% 98%
2
53% 73%
56% 32% 75%
19%
88%
4
Number of nodes
8
16
(c) Weak Scaling 95%
93% 95%
93% 100%
2 N=4 (LX2) N=7 (LX2)
94% 95%
99%
4
8
N=15 (LX2) N=4 (A100)
N=7 (A100)
Number of nodes
16
Fig. 10. Intra/Inter-node Scaling Experiments
of excessive shared-memory demand, whereas SME reaches near-full utilization and delivers approximately 1.6× higher FLOP rate than at lower order. Our cross-platform results show that SME materially improves the competitiveness of ARM multicore CPUs for tensor-product SEM, reducing the gap to accelerator-oriented platforms and changing the practical role of CPUs from baseline hosts to viable highperformance execution targets. Graph coloring reduces the number of elements processed within each parallel region, making performance more sensitive to working-set size; if the local workset becomes too small, hardware resources cannot be fully saturated. From this perspective, the higher-order cases in Sec.V-C have not yet reached their full performance potential, and the gain from larger p would likely increase further for larger worksets. For the inter-node weak-scaling study, we report only p = 7 and p = 15, since the p = 4 case requires too many elements and exceeds the capability of our mesher. The weak-scaling results are close to ideal. As discussed for the communication-
thread optimization, the halo region is small under the hybrid execution model: each process exchanges less than 1 MB of data compared with a working set of several gigabytes. As a result, communication overhead is minor, and the observed deviation from ideal weak scaling is dominated instead by system noise and by the 5–10% load imbalance introduced by the METIS/Scotch partitioning. For strong scaling, the problem size is fixed while the number of resources is increased, so the work per process decreases and parallel overhead becomes more visible. Under this regime, p = 15 exhibits the weakest strong scaling because our optimizations are sensitive to local element count and working-set size, whereas p = 4 shows the strongest strong-scaling behavior. This is consistent with the standard observation that strong scaling degrades once the local workload becomes too small to keep the hardware fully utilized. VI. C ONCLUSION In this work, we optimized SPECFEM3D for ARM multicore CPUs with SME and showed that the benefit of SME extends beyond faster tensor-product kernels. By combining an SME-aware batched small-matrix implementation with a memory-aware hybrid MPI+OpenMP execution scheme, we substantially improved full-application performance on LX2 while enabling realistic workloads under limited HBM capacity. More importantly, through a dispersion-based isoaccuracy study, we showed that SME changes the practical cost trend across polynomial orders and makes higher-p SEM discretizations significantly more attractive on modern ARM multicore CPUs. These results suggest that matrix-oriented CPU extensions can affect not only kernel mapping, but also application-level discretization choices in high-order wavepropagation codes. R EFERENCES [1] Z. Zhang, W. Zhang, and X. Chen, “Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics,” Geophysical Journal International, vol. 199, no. 2, pp. 860–879, 2014. [2] B. Chen, H. Fu, Y. Wei, C. He, W. Zhang, Y. Li, W. Wan, W. Zhang, L. Gan, W. Zhang, Z. Zhang, G. Yang, and X. Chen, “Simulating the wenchuan earthquake with accurate surface topography on sunway taihulight,” ser. SC ’18. IEEE Press, 2018. [3] W. Wan, L. Gan, W. Wang, Z. Yin, H. Tian, Z. Zhang, Y. Wang, M. Hua, X. Liu, S. Xiang, Z. He, Z. Wang, P. Gao, X. Duan, W. Liu, W. Xue, H. Fu, G. Yang, X. Chen, Z. Song, Y. Chen, X. Liu, and W. Zhang, “69.7-pflops extreme scale earthquake simulation with crossing multi-faults and topography on sunway,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’23. New York, NY, USA: Association for Computing Machinery, 2023. [Online]. Available: https://doi.org/10.1145/3581784.3613209 [4] H.-W. Zhou, H. Hu, Z. Zou, Y. Wo, and O. Youn, “Reverse time migration: A prospect of seismic imaging methodology,” EarthScience Reviews, vol. 179, pp. 207–227, 2018. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0012825217306256 [5] C. Li, G. Liu, and Y. Li, “A practical implementation of 3d tti reverse time migration with multi-gpus,” Comput. Geosci., vol. 102, no. C, p. 68–78, May 2017. [Online]. Available: https://doi.org/10.1016/j.cageo.2017.02.011 [6] D. Komatitsch and J. Tromp, “Introduction to the spectral element method for three-dimensional seismic wave propagation,” Geophysical Journal International, vol. 139, no. 3, pp. 806–822, 12 1999. [Online]. Available: https://doi.org/10.1046/j.1365-246x.1999.00967.x
[7] ——, “Spectral-element simulations of global seismic wave propagation—i. validation,” Geophysical Journal International, vol. 149, no. 2, pp. 390–412, 05 2002. [Online]. Available: https://doi.org/10.1046/j.1365-246X.2002.01653.x [8] C. G. Canuto, M. Y. Hussaini, A. M. Quarteroni, and T. A. Zang, Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics (Scientific Computation). Berlin, Heidelberg: SpringerVerlag, 2007. [9] T. Germann, T. Kolev, P. Fischer, M. Min, J. Dongarra, J. Brown, V. Dobrev, T. Warburton, S. Tomov, M. S. Shephard, A. Abdelfattah, V. Barra, N. Beams, J.-S. Camier, N. Chalmers, Y. Dudouit, A. Karakus, I. Karlin, S. Kerkemeier, Y.-H. Lan, D. Medina, E. Merzari, A. Obabko, W. Pazner, T. Rathnayake, C. W. Smith, L. Spies, K. Swirydowicz, J. Thompson, A. Tomboulides, and V. Tomov, “Efficient exascale discretizations: High-order finite element methods,” vol. 35, no. 6, p. 527–552, Nov. 2021. [Online]. Available: https://doi.org/10.1177/10943420211020803 [10] K. Swirydowicz, N. Chalmers, A. Karakus, and T. Warburton, “Acceleration of tensor-product operations for high-order finite element methods,” Int. J. High Perform. Comput. Appl., vol. 33, no. 4, p. 735–757, Jul. 2019. [Online]. Available: https://doi.org/10.1177/1094342018816368 [11] D. Peter, D. Komatitsch, Y. Luo, R. Martin, N. Le Goff, E. Casarotti, P. Le Loher, F. Magnoni, Q. Liu, C. Blitz, T. Nissen-Meyer, P. Basini, and J. Tromp, “Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes,” Geophysical Journal International, vol. 186, no. 2, pp. 721–739, 08 2011. [Online]. Available: https://doi.org/10.1111/j.1365-246X.2011.05044.x [12] J. Andrej, N. Atallah, J.-P. Bäcker, J.-S. Camier, D. Copeland, V. Dobrev, Y. Dudouit, T. Duswald, B. Keith, D. Kim, T. Kolev, B. Lazarov, K. Mittal, W. Pazner, S. Petrides, S. Shiraiwa, M. Stowell, and V. Tomov, “High-performance finite elements with mfem,” The International Journal of High Performance Computing Applications, Jun. 2024. [Online]. Available: http://dx.doi.org/10.1177/10943420241261981 [13] S. Touhami, F. Gatti, F. Lopez-Caballero, D. A. C. Cruz, and D. Clouteau, “SEM3D: A 3d high-fidelity numerical earthquake simulator for broadband (0–10 hz) seismic response prediction at a regional scale,” Geosciences, vol. 12, no. 3, p. 112, 2022. [14] G. Bell, D. H. Bailey, J. Dongarra, A. H. Karp, and K. Walsh, “A look back on 30 years of the gordon bell prize,” Int. J. High Perform. Comput. Appl., vol. 31, no. 6, p. 469–484, Nov. 2017. [Online]. Available: https://doi.org/10.1177/1094342017738610 [15] S. Henneking, S. Venkat, V. Dobrev, J. Camier, T. Kolev, M. Fernando, A.-A. Gabriel, and O. Ghattas, “Real-time bayesian inference at extreme scale: A digital twin for tsunami early warning applied to the cascadia subduction zone,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’25. New York, NY, USA: Association for Computing Machinery, 2025, p. 60–71. [Online]. Available: https://doi.org/10.1145/3712285.3771787 [16] S. Matsuoka, M. Sato, H. Kobayashi et al., “Co-design and system for the supercomputer “fugaku”,” IEEE Micro, 2022. [17] RIKEN, “Riken launches international initiative with fujitsu and nvidia to accelerate development of fugakunext,” https://www.riken.jp/en/news pubs/news/2025/20250822 1/index.html, 2025, accessed: 2026-04-07. [18] Arm Ltd., “Arm architecture reference manual for a-profile architecture,” https://developer.arm.com/documentation/ddi0487/latest/, 2024, includes the Scalable Matrix Extension (SME); accessed 2026-04-07. [19] D. Arndt, W. Bangerth, D. Davydov, T. Heister, L. Heltai, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin, and D. Wells, “The deal.ii finite element library: Design, features, and insights,” Computers & amp; Mathematics with Applications, vol. 81, p. 407–422, Jan. 2021. [Online]. Available: http://dx.doi.org/10.1016/j.camwa.2020.02.022 [20] A. Abdelfattah, V. Barra, N. Beams, R. Bleile, J. Brown, J.-S. Camier, R. Carson, N. Chalmers, V. Dobrev, Y. Dudouit, P. Fischer, A. Karakus, S. Kerkemeier, T. Kolev, Y.-H. Lan, E. Merzari, M. Min, M. Phillips, T. Rathnayake, R. Rieben, T. Stitt, A. Tomboulides, S. Tomov, V. Tomov, A. Vargas, T. Warburton, and K. Weiss, “Gpu algorithms for efficient exascale discretizations,” 2021. [Online]. Available: https://arxiv.org/abs/2109.05072 [21] A. Heinecke, G. Henry, M. Hutchinson, and H. Pabst, “Libxsmm: accelerating small matrix multiplications by runtime code generation,” in Proceedings of the International Conference for High Performance
Computing, Networking, Storage and Analysis, ser. SC ’16. IEEE Press, 2016. [22] M. T. Jones and P. E. Plassmann, “A parallel graph coloring heuristic,” SIAM Journal on Scientific Computing, vol. 14, no. 3, pp. 654–669, 1993. [Online]. Available: https://doi.org/10.1137/0914041 [23] P. Krysl, “Parallel assembly of finite element matrices on multicore computers,” Computer Methods in Applied Mechanics and Engineering, vol. 428, p. 117076, 2024. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0045782524003323 [24] C. Cecka, A. Lew, and E. Darve, “Assembly of finite element methods on graphics processors,” International Journal for Numerical Methods in Engineering, vol. 85, pp. 640 – 669, 02 2011. [25] D. Komatitsch, G. Erlebacher, D. Göddeke, and D. Michéa, “High-order finite-element seismic wave propagation modeling with mpi on a large gpu cluster,” Journal of Computational Physics, vol. 229, pp. 7692– 7714, 10 2010. [26] R. Rabenseifner, G. Hager, and G. Jost, “Hybrid mpi/openmp parallel programming on clusters of multi-core smp nodes,” in 2009 17th Euromicro International Conference on Parallel, Distributed and Networkbased Processing, 2009, pp. 427–436. [27] G. Schubert, H. Fehske, G. Hager, and G. Wellein, “Hybrid-parallel sparse matrix-vector multiplication with explicit communication overlap on current multicore-based systems,” CoRR, vol. abs/1106.5908, 2011. [Online]. Available: http://arxiv.org/abs/1106.5908 [28] E. Castillo, N. Jain, M. Casas, M. Moreto, M. Schulz, R. Beivide, M. Valero, and A. Bhatele, “Optimizing computation-communication overlap in asynchronous task-based programs,” in Proceedings of the ACM International Conference on Supercomputing, ser. ICS ’19. New York, NY, USA: Association for Computing Machinery, 2019, p. 380–391. [Online]. Available: https://doi.org/10.1145/3330345.3330379