arXiv CS open access

Beyond Fast Contractions: Attenuation and Recovery of Matrix-Engine Speedups in High-Order Finite Elements

2026 · arxiv_cs

arXiv CS · Papers · 2026

cloud distributed-computing parallel-computing

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

arXiv:2608.09175v1 [cs.DC] 10 Aug 2026

Beyond Fast Contractions: Attenuation and Recovery of Matrix-Engine Speedups in High-Order Finite Elements Yinuo Wang

Lin Gan

Tianqi Mao

wyn22@mails.tsinghua.edu.cn Tsinghua University Beijing, China

lingan@tsinghua.edu.cn Tsinghua University Beijing, China Hetao Institute of Mathematics and Interdisciplinary Sciences Shenzhen, China

mtq24@mails.tsinghua.edu.cn Tsinghua University Beijing, China

Zeyu Song

Wubing Wan

Jiayu Fu

zy-song22@mails.tsinghua.edu.cn Tsinghua University Beijing, China

wanwubing@himis-sz.cn Hetao Institute of Mathematics and Interdisciplinary Sciences Shenzhen, China

fujy24@mails.tsinghua.edu.cn Tsinghua University Beijing, China

Zekun Yin

Yuyang Jin

Xiaohui Duan

Zekun.yin@sdu.edu.cn Shandong University Jinan, China

jinyuyang@tsinghua.edu.cn Tsinghua University Beijing, China

sunrise.duan@sdu.edu.cn Shandong University Jinan, China

Wei Xue

Guangwen Yang

xuewei@tsinghua.edu.cn Tsinghua University Beijing, China

ygw@tsinghua.edu.cn Tsinghua University Beijing, China

Abstract

Modern processors increasingly provide matrix engines whose peak arithmetic throughput greatly exceeds conventional SIMD, but scientific applications rarely realize this advantage end to end. We examine this gap in SPECFEM3D’s dominant stiffness operator on the Arm LX2 CPUs that power the flagship Lineshine supercomputer. Against a matched, high-performance SVE baseline on the same cores, SME’s 4× single-precision peak advantage falls to 2.2× for isolated tensor contractions and 1.1× for the complete operator. Our factorized diagnostic attributes the loss to pointwise computation, indirect field movement and synchronization, and irregular coefficient delivery. Explicit SIMD mitigates pointwise work, raising the full-operator speedup to 1.3×. Field-layout changes mitigate indirect movement and synchronization, while vector-blocked coefficient streaming reduces irregular-access costs; together they raise speedup to 1.6× at high order. A contraction-free control bounds further contraction-only gains at 1.11–1.32×. Realizing matrix-engine performance therefore requires co-designing the entire operator path, not merely replacing its contraction kernel.

On modern manycore and accelerator-based systems, data movement is often substantially more expensive than arithmetic in both time and energy. Matrix-free high-order finite and spectral element methods (FEM/SEM) avoid assembled sparse matrices. They exploit tensor-product elements, replacing assembled-matrix accesses with local computation. This reduces global memory traffic while exposing high arithmetic intensity and fine-grained parallelism. These properties improve accuracy per computational time and make high-order discretizations important for modern manycore and accelerator-based systems [1, 23], including wave propagation, high-order computational fluid dynamics (CFD), and atmospheric modeling [12, 24, 25, 29, 34]. On the hardware side, matrix units are a major source of peak floating-point capability, made prominent by GPU Tensor Cores for dense GEMM. Recent work remaps sparse products and broader scientific patterns onto them [15, 27, 31, 40]. In high-order stiffness, sum factorization expresses tensor-product differentiation and weak-form application as batches of small contractions, making contraction a natural matrix-unit target [11, 35]. This apparent hardware–algorithm match is incomplete. A highorder stiffness operator is not a single contraction kernel: pointwise physics and metrics, field gathers and scatters, reconciliation of shared degrees of freedom (DoFs), and irregular coefficient access surround the contractions. These stages can attenuate matrix-unit

Keywords Arm SME, spectral element method, high-order finite elements, matrix-free operators, tensor contractions, data layout, performance analysis 1

Introduction

107

(a) Lineshine weak scaling Forward iteration time (%)

Estimated irregular-elastic GFLOP/s

Yinuo Wang, Lin Gan, Tianqi Mao, Zeyu Song, Wubing Wan, Jiayu Fu, Zekun Yin, Yuyang Jin, Xiaohui Duan, Wei Xue, and Guangwen Yang

p=4 p=7 p = 15

106 105 104 20

24 26 Nodes

210

100

(b) Runtime composition

This cascade is the starting point of the paper. Accordingly, we ask: relative to a matched, high-performance single-instruction, multiple-data (SIMD) baseline, where is matrix-unit acceleration exposed or masked across a complete high-order SEM stiffness operator, and which operator-path and data-layout changes recover it? To study this cascade, we propose a factorized diagnostic methodology for high-order FEM stiffness operators that separates operatorlevel decomposition from architecture-specific execution backends. It independently controls the contraction backend, pointwise physics and geometry implementation, and field and coefficient layouts while comparing matched null, SIMD, and matrix-unit contraction backends. Comparisons among these controlled variants use measured time, throughput, and effective-bandwidth deltas—rather than architecture-specific hardware counters—to identify where matrix-unit acceleration is exposed, masked, or lost. We instantiate the methodology on LX2 using SVE and SME and apply it to SPECFEM3D’s acoustic and elastic formulations [30]. Elastic is our primary coupled vector-field target, while acoustic provides a scalar control for understanding general FEM stiffness behavior.2 Porting the methodology to another matrix-engine CPU preserves the same decomposition and experimental controls while replacing the SVE and SME implementations with native vector and matrix-unit backends. The first target is the conventional SPECFEM3D-style L-vector operator. This shared-DoF layout is compact and widely used, but surrounds each contraction with indirect gathers, shared accumulations, and race-avoidance constraints. The diagnosis first identifies a fixable attenuation layer: pointwise physics and metric work can hide backend differences unless controlled with explicit SIMD. After this fix, attenuation remains in the inherited L-vector structure: shared DoFs create indirect access and race-avoidance costs, and irregular geometry adds coefficient and metric traffic that must be separated from contraction time. We then use the diagnosis to guide mitigation. For field movement, S-vector and E-vector layouts move shared-DoF costs to different parts of the execution path: S-vector reduces conflicts through subdomain duplication, while E-vector exposes regular element-local kernels at the cost of explicit reconciliation. The choice is a tradeoff between indirect access and scatter conflicts inside the kernel, or explicit synchronization after a more regular local kernel. For irregular coefficients, array-of-structures-of-vectors (AoSoV) and mixed-AoSoV layouts make coefficient access more vector-friendly and reduce traffic. In the elastic irregular case, this recovery path raises the high-order SME/SVE operator speedup from about 1.1× in the original L-vector diagnostic to about 1.6× after field and coefficient layout recovery. The stronger recovered gains at high polynomial order motivate a smooth, high-fidelity wave-propagation case study of how architecture-dependent operator throughput can influence equal-accuracy discretization choices. Taken together, the attenuation diagnosis and layout-mediated recovery reveal the paper’s central insight: matrix-engine acceleration is governed not by the contraction backend alone, but by where

80 60 40 20 0

Stiffness Comm.

Mass inverse Newmark update

p=4 p=7 p = 15 512 nodes 1024 nodes 1024 nodes

Figure 1: SPECFEM3D weak scaling and solver runtime composition on Lineshine for polynomial orders 𝑝 ∈ {4, 7, 15}. The application retains high parallel efficiency while stiffness remains the dominant measured hotspot at scale.

LX2 Peak FP32 Ratio(4x)

Hardware

Contraction

Operator

SME-SVE parity(1x)

Figure 2: Peak-to-operator attenuation of SME over SVE. Across polynomial order, horizontal lines mark the LX2 FP32 matrix-FMA peak ratio (4×) and parity (1×); the measured curves expose the successive loss from matched contractions to the conventional compiler-vectorized irregular-elastic Lvector operator.

gains even when the contractions themselves map efficiently, creating a gap between matrix-unit peak, contraction performance, and complete-operator speedup. This gap becomes consequential as matrix units move into generalpurpose multicore CPUs. Arm’s Scalable Matrix Extension (SME) is an important example. Lineshine pushes this trend to system scale, sustaining 2.198 EFLOP/s on HPL as the top-ranked June 2026 TOP500 system [33] with LX2 CPUs whose single-precision peak is dominated by SME [6]. Lineshine provides a scalable production case: SPECFEM3D retains about 94–98% weak-scaling efficiency through 512–1024 nodes, with stiffness consuming 59–79% of forward-iteration time at the largest scales (Figure 1). Yet SME’s 4× FP32 matrix-FMA peak advantage over Arm’s Scalable Vector Extension (SVE) falls to 2.2× for isolated contractions against a matched high-performance SVE proxy. The full irregular-elastic L-vector1 operator reaches only 1.1× (Figure 7). We call this loss from matrix-unit peak through contractions to full stiffness the peak-to-operator attenuation cascade (Figure 2).

2 For space, detailed figures focus on elastic, while the roofline retains both physics.

Complete measurements, benchmark source code, data, and plotting scripts are available at https://github.com/FrankTrek/SME_HighOrderFEM_Benchmark.

1We use libCEED terminology [8]; see Sections 2 and 6.1 for details.

Beyond Fast Contractions: Attenuation and Recovery of Matrix-Engine Speedups in High-Order Finite Elements Die 1

...

Continuous-Galerkin Elements (Rank-Local) Element

Vector B

Vector A

Die 0

ZA Tile Cores

HBM

DDRC

SME Matrix Unit

LX2 CPU

Figure 3: LX2 organization and per-core SME execution. The dual-die package contains four 38-core NUMA domains per die, connected to HBM stacks, DDR controllers, and I/O. Every CPU core contains an SME matrix unit. The inset illustrates the outer product of vectors A and B accumulated into ZA using 𝑉𝐿 = 8 for visual clarity (our measured FP32 configuration uses 𝑉𝐿 = 16).

55 56

Mapping from Element-Local Node to Shared Vector Entry Element Restriction

Indirect Gather

Element-Local Tensor

Ordinary (non-shared) nodes Shared nodes between elements

Figure 4: Conventional L-vector storage and element restriction. A rank-local continuous-Galerkin vector assigns one entry to each distinct DoF, with interface entries shared by neighboring elements. An indirect restriction map gathers these entries into a dense element-local tensor for operator evaluation.

the complete operator pays for pointwise work, field movement, coefficient delivery, and shared-DoF reconciliation. This paper makes four contributions: (1) Matched SVE/SME contraction backends with an attainable performance bound. We build strong vectorunit and matrix-unit tensor-contraction backends and show that SME approaches the order-dependent bound for small SEM contractions. (2) A factorized full-operator diagnosis of matrix-unit attenuation. Matched null, SVE, and SME backends and controlled operator variants expose how pointwise physics, field movement, coefficient delivery, and reconciliation attenuate acceleration from contraction to the complete operator. (3) Layout-mediated recovery with an operator-level envelope. Explicit SIMD, S/E-vector, and AoSoV variants recover successive attenuation layers, reaching about 1.6× high-order elastic-irregular SME/SVE speedup; the roofline spans both physics and geometries. (4) A discretization implication for high-order FEM. We find evidence that the recovered high-order operator throughput favors higher polynomial order in the tested equalaccuracy wave-propagation setting.

outer product is accumulated into the architectural ZA tile, and data moves between ZA and SVE registers through horizontal or vertical slice instructions. ZA consists of 64 rows of 64 bytes, forming a twodimensional architectural tile, so practical SEM contractions occupy only part of the tile at low polynomial order. For single-precision FMA, SME provides approximately 4× the theoretical peak throughput of conventional SVE multiply-accumulate instructions on the same core. This peak ratio motivates the comparison, but the rest of the paper studies how much survives once the contraction is embedded in a full operator.

2.2

High-Order FEM and Tensor Contraction

The stiffness operators arise from weak forms of acoustic and elastic wave propagation. For scalar acoustic propagation, we use ∫ 𝑎 ac (𝑢, 𝑣) = 𝜌 −1 ∇𝑣 · ∇𝑢 d𝑥, (1) Ω

where 𝑢 and 𝑣 are scalar trial and test fields and 𝜌 −1 is inverse density. For elastic propagation, we use ∫ 𝑎 el (u, v) = ∇v : C : ∇u d𝑥, (2)

2 Background 2.1 The LX2 CPU Architecture

where u and v are displacement trial and test fields and C is the elastic stiffness tensor. In the isotropic case used here, C is parameterized by the Lamé parameters 𝜆 and 𝜇. Each physical hexahedral element Ω𝑒 is represented by a map b = [−1, 1] 3 ; its x𝑒 (𝜉, 𝜂, 𝜁 ) = (𝑥, 𝑦, 𝑧) from the reference cube Ω tensor-product Gauss–Lobatto–Legendre (GLL) nodes lie at (𝜉𝑖 , 𝜂 𝑗 , 𝜁𝑘 ). SEM evaluates gradients through this tensor-product structure. At polynomial order 𝑝, let 𝑛 = 𝑝 + 1 denote the number of GLL points in each reference-coordinate direction. The element-local field is a small three-dimensional tensor 𝑢𝑖 𝑗𝑘 , with 𝑖, 𝑗, 𝑘 = 0, . . . , 𝑛 − 1. Let 𝐷 denote the one-dimensional differentiation matrix evaluated at the GLL points. The three reference-space derivatives differ only in the

The LX2 processor used in Lineshine is a dual-die Armv9-A CPU with 304 physical cores running at 1.6 GHz (Figure 3). Each die is divided into four NUMA domains of 38 cores. Every core provides conventional SVE/SVE2 vector execution and supports SME, while the memory system combines per-NUMA high-bandwidth memory (HBM) with socket-level DDR. Each NUMA domain provides 4 GB of HBM, which reaches about 400 GB/s in our STREAM measurement [28], while DDR provides larger capacity at about 120 GB/s. SME evaluates matrix-style computation through an outer-product execution model. Two vectors are read from SVE registers, their 3

Yinuo Wang, Lin Gan, Tianqi Mao, Zeyu Song, Wubing Wan, Jiayu Fu, Zekun Yin, Yuyang Jin, Xiaohui Duan, Wei Xue, and Guangwen Yang

tensor dimension to which 𝐷 is applied; using 𝜉 as representative, (𝜕𝜉 𝑢)𝑖 𝑗𝑘 =

𝑛−1 ∑︁

𝐷𝑖𝑎𝑢𝑎 𝑗𝑘 .

(3)

𝑎=0

The 𝜂 and 𝜁 forms apply 𝐷 to the second and third indices, respectively. Each directional derivative is a tensor contraction, equivalently a batch of small matrix multiplications over element-local slices [9]. The test-function gradients enter through a transposegradient operation with the same tensor-product structure using 𝐷𝑇 . Between the gradient and transpose-gradient contractions, the element kernel performs pointwise physics and geometry work. At each quadrature point, a 3 × 3 contravariant matrix maps referencespace derivatives to physical space. The acoustic factor 𝜌 −1 or the elastic tensor C is applied. The resulting physical flux is mapped back to reference space before the transpose-gradient forms the element-local residual. The L-vector maps each element-local GLL node to a shared vector entry (Figure 4). Its compact storage requires indirect gather/scatter, while concurrent updates require atomics or conflict-free element coloring [21]. Restriction, directional contractions, pointwise physics and geometry, transpose contractions, and accumulation define the attenuation path studied below. Section 4 calibrates the matrix-unit opportunity in the contractions. Section 5 separates pointwise physics, field movement, and coefficient traffic in the conventional L-vector operator. Section 6 changes field and coefficient representations to recover the exposed losses.

Experimental Setup, Validation, and Data Policy

All measurements use the vendor-provided LX2 platform. Each run is confined to one HBM NUMA domain and uses 36 CPU cores, with HBM holding the hot stiffness working set. OpenMP threads are pinned to physical cores with OMP_PLACES=cores and OMP_PROC_BIND=close; each timed invocation follows 50 warmup iterations with 200 measured iterations. We report FP32 unless otherwise stated. SVE and SME binaries are built from matched kernel contracts, and every optimized path is validated through canonical operator outputs against the native flat-L-vector reference (10−4 absolute and relative tolerances). Mixed-AoSoV storage is accepted only when relative 𝐿2 and 𝐿∞ errors are each at most 10−3 . Builds use clang++ with -mcpu=hip11, -O3, -ffast-math, and a fixed 512-bit SVE vector length. The explicit-SIMD design-space sweep provides 16 process-level samples for each of 840 configurations, yielding 13,440 timings; the compiler-vectorized control uses five. A scripted stability check identifies 246 samples (1.83%) as intermittent high-side slowdowns. These samples are sparse in most configurations and are excluded before taking the median. For the few problematic configurations, the slowdowns affect a substantial fraction of repetitions or exhibit collection-cohort dependence, making filtering alone unreliable. We therefore remeasure these configurations using unchanged binaries and matched null, SVE, and SME commands executed sequentially on one dedicated HBM NUMA domain. We report the median of 4

five complete process-level repetitions without further filtering, retaining a point only when all outputs validate. The artifact repository contains original and rerun measurements, command lines, binary hashes, sweep manifests, filtering code, validation tolerances, NUMA placement, OpenMP settings, compiler flags, thresholds, exclusions, and backend build macros. The manuscript reports the selected summary points needed for the performance narrative.

Contraction-Layer Attenuation: Bound and Matched Backends

The first attenuation layer appears before the full stiffness operator is considered. Both backends evaluate the directional gradient and transpose-gradient contractions represented by Equation (3): SVE uses unrolled, register-tiled SIMD FMA kernels, while SME maps the work to outer-product updates on ZA tiles. Despite SME’s higher peak, SEM contractions are small. Each directional derivative acts on length-𝑛 lines inside an 𝑛 × 𝑛 × 𝑛 cube, so each line exposes only 𝑛 useful contraction opportunities. Before 𝑛 approaches the SME vector/tile width 𝑉𝐿 , the attainable utilization of matrix-FMA peak scales as 𝑛/𝑉𝐿 . This is a shape bound imposed by the operator, not a code-generation artifact. Load/store overhead further lowers the attainable rate. For the 𝜂- and 𝜁 -direction contractions, data can be aggregated into SMEfriendly operands, but accumulated results still have to be stored. We model the first-order bound as 𝑛 CPI OPA · FLOPS SME . (4) CPI OPA + CPI LSU 𝑉𝐿 Here CPI denotes cycles per instruction; OPA and LSU denote outerproduct accumulation and the load/store unit, respectively. The 𝜉 direction has an additional layout penalty because element-local memory order does not match the operand layout required by outer products. The kernel performs in-register transpose and reorganization to avoid strided loads. This adds roughly one extra load and store stream, lowering the bound to CPI OPA 𝑛 · FLOPS SME . (5) CPI OPA + 3CPI LSU 𝑉𝐿 For the 512-bit FP32 configuration, 𝑉𝐿 = 16 lanes and the rounded single-core matrix-FMA peak is FLOPS SME = 400 GFLOP/s. The modeled issue costs CPI OPA = 2 and CPI LSU = 1 yield attainablerate factors of 2/3 for 𝜂, 𝜁 and 2/5 for 𝜉. Our backend is designed to approach these direction-specific bounds (Figure 5). It aggregates work across the batched slice dimension and pipelines independent slices across ZA tiles. For the 𝜁 direction, element layout exposes contiguous 𝜉𝜂 planes, so the kernel aggregates them into a larger 𝑛 × 𝑛 2 multiplication. The 𝜉 direction uses in-register transpose/reorganization to recover contiguous vector access. The 𝜂 direction is closest to a direct batched small matrix multiplication over 𝜁 slices; for 𝑝 = 4–7, it additionally uses vector concatenation and shuffle operations to pack multiple short element-local lines into one outer-product operand. The matched SVE path ensures that later measurements compare hardware datapaths rather than a strong SME implementation against a weak SIMD baseline. Figure 6 reports single-core microbenchmarks. At high order, the SVE backend reaches about 65 GFLOP/s, roughly 65% of the

Beyond Fast Contractions: Attenuation and Recovery of Matrix-Engine Speedups in High-Order Finite Elements

η-Direction Derivative

ξ-Direction Derivative

Data Layout in Memory

Nodal Values Slices

Differentiation Matrix

For N ≤ VL/2, one SIMD register can hold two η-lines

Load Vector η-th Line

Load Vector

(η+1)-th Line

(ζ, η)

(ζ, η+1)

(ζ+1, η)

(ζ+1, η+1)

Horizontal Load

(ζ+1)-th Slice

TBL shuffle

ζ-th Slice

Derivatives

In Memory

(ζ, η)

(ζ+1, η)

(ζ, η+1)

(ζ+1, η+1)

ZA Tile

Transpose Required Vector Discontinuous

ζ-Direction Derivative Load Vector Extract Slice along η-axis

In Memory

Naturally Concatenated (ζ, η)

(ζ, η+1)

Data Layout in Memory

Figure 5: Direction-specific data movement in the SME contraction backend. The 𝜂 path packs short lines, the 𝜉 path reorganizes discontinuous memory operands through ZA-tile transposition, and the 𝜁 path exploits naturally contiguous slices.

(a) ξ-summation Performance (GFLOP/s)

250

(b) η-summation

Theoretical CPI Bound

(d) Full contraction kernel

200 150

Theoretical CPI Bound

100 50 0

4 5 6 7

Polynomial order p

4 5 6 7

Polynomial order p Compiler-Opt

SVE

4 5 6 7

Polynomial order p

SME

4 5 6 7

Polynomial order p

Directional Harmonic Mean

Figure 6: Single-core directional tensor contractions and full-contraction proxy. The SME curves should be compared with the order- and direction-dependent attainable bounds rather than nominal matrix-FMA peak. single-core SVE FMA peak; given the small matrix sizes and pipeline and load/store overheads, it is a strong vector-engine baseline. The directional results show that the CPI-style bound is the right point of comparison. Performance grows as 𝑛 occupies more of ZA, while the 𝜂 and 𝜁 curves retain visible store overhead and 𝜉 remains lower because of transpose/reorganization. The 𝜂 direction drops after 𝑝 = 7, where the low-order concatenation path no longer applies, but resumes the same tile-utilization trend at higher orders. Across directional kernels, SME reaches up to 250 GFLOP/s and outperforms SVE by up to 3.6× at high order. At 𝑝 = 15, the speedups are 2.1×, 3.6×, and 3.3× for 𝜉, 𝜂, and 𝜁 , respectively. Compared with the compiler-optimized scalar backend, the directional SME kernels are one order of magnitude faster at high order.

The full contraction proxy is the relevant calibration for stiffness. Because all three directions execute in sequence, its effective rate is governed by the slower direction through a harmonic-mean effect. At 𝑝 = 15, the full-proxy SME/SVE speedup is about 2.2×, close to the 𝜉 speedup and well below 𝜂 and 𝜁 . The compiler-optimized backend is not a reliable comparison target: it changes behavior sharply across orders and exhibits a cliff at 𝑝 = 15. The gap between the directional harmonic mean and full proxy also grows at high order, consistent with the working set exceeding the 32 KB L1 cache. Thus a 4× hardware peak advantage becomes about 2.2× at the contraction layer even before full-operator costs are introduced.

Yinuo Wang, Lin Gan, Tianqi Mao, Zeyu Song, Wubing Wan, Jiayu Fu, Zekun Yin, Yuyang Jin, Xiaohui Duan, Wei Xue, and Guangwen Yang

L-vector Attenuation Layers

Explicit-SIMD pointwise path

SVE

SME

Pointwise Physics Auto-Vec

Intrinsic

Field Movement

Coefficient Traffic

Irregular SoA coefficient traffic (a) Compiler-vectorized, regular

Contraction Backend Null

Controlled L-Vector Attenuation Sweep

(b) Compiler-vectorized, irregular SoA

(d) Explicit SIMD, irregular SoA

Contraction–Memory Overlap Diagnostic

Effective Bandwidth Delivery

p=15 p=13 p=11 p=9 p=7 p=5 p=4

0.5

1.0

Time normalized to null irregular

1.5

Figure 7: Factorized diagnosis of L-vector attenuation for elastic stiffness. The sidebar traces contraction backend, pointwise SIMD, field movement, and irregular SoA coefficient traffic. (A) Controlled backend sweeps separate SIMD and geometry effects. (B) The partial-overlap analysis compares observed, null-backend, and estimated contraction costs. (C) Effective-bandwidth measurements expose the residual shared-DoF and coefficient-streaming limitations. (3) Total null-backend bandwidth. As a conservative aggretotal = (𝑄 gate, we report 𝐵 null var + 𝑄 coeff )/𝑇null . To reason about why SVE or SME does not simply add contrac𝑏 tion time to the null time, we use a partial-overlap model. Let 𝑇bcontr be the contraction time estimated from Section 4. We write

5 L-Vector Attenuation Diagnosis 5.1 Diagnostic Design and Metrics The benchmark preserves the matrix-free operator while omitting file, mesh-I/O, and orchestration infrastructure so implementations can vary independently. Its axes are contraction backend (null, SVE, SME), regular or irregular geometry, scalar acoustic or vector elastic physics, and compiler-vectorized or explicit-SIMD pointwise work. Together they separate contraction, variable movement, physics, and coefficient traffic. Let 𝑇null denote null-backend stiffness time and 𝑇𝑏 the measured time for 𝑏 ∈ {SVE, SME}. The modeled hot-path traffic is split into variable-field traffic 𝑄 var and coefficient/metric traffic 𝑄 coeff . We use three diagnostic quantities:

𝑏 𝑏 𝑇𝑏 ≈ max(𝑇mem, 𝑇bcontr ) + 𝛼𝑏 min(𝑇mem, 𝑇bcontr ),

(6)

where 𝑇mem is null-backend memory/control time and 𝛼𝑏 is a fitted non-overlap parameter. Values near zero indicate good overlap; values near one indicate near serialization. The model is diagnostic rather than cycle accurate. Table 1 uses 2.02–2.59 GiB elastic irregular-SoA working sets, chosen to fit within one NUMA domain’s HBM while providing enough independent element work to occupy all 36 cores and measure steady-state, resource-saturated performance. From 𝑝 = 4 to 15, contraction length grows from 𝑛 = 5 to 16 and modeled FLOP share from 60.4% to 83.0%, marking the transition from shortline, larger-pointwise-share cases to full-width, contraction-heavier cases. These shares are not runtime fractions.

(1) Variable-field bandwidth. For regular geometry, the null reg backend gives the control 𝐵 var = 𝑄 var /𝑇null . (2) Incremental coefficient bandwidth. For irregular geometry, the irregular-minus-regular delta estimates pointwiseextra = 𝑄 irreg /(𝑇 irreg − 𝑇 reg ). This is an stream cost as 𝐵 coeff coeff null null effective incremental quantity, not raw copy bandwidth: the timing subtraction can include overlap with variable movement and arithmetic.

5.2

Attenuation in the L-Vector Path

We apply the L-vector attenuation diagnosis to both physics. For elastic, Figure 7 shows the complete factorized diagnosis for our 6

Beyond Fast Contractions: Attenuation and Recovery of Matrix-Engine Speedups in High-Order Finite Elements

Table 1: Elastic irregular-SoA L-vector workload. All cases use ℎ = 0.5; mesh dimensions vary to keep the HBM-resident problem scale comparable. Footprint includes field, numbering, coefficient, plan, reference-element, and scratch storage. Work is analytical GFLOP per application, and contraction share is a modeled FLOP fraction rather than measured time. 𝑝

Mesh

𝑁𝑒 (k)

MDoF

GiB

GF/app.

Contr. (%)

4 5 7 9 11 13 15

63×63×56 54×54×44 42×42×36 36×33×28 27×27×24 24×24×20 24×21×16

222.3 128.3 63.5 33.3 17.5 11.5 8.1

43.21 48.69 66.05 73.51 70.60 76.71 82.48

2.02 2.08 2.43 2.55 2.33 2.47 2.59

8.28 9.26 13.20 15.90 16.63 19.66 22.92

60.4 64.7 70.9 75.3 78.5 81.0 83.0

locality sequence as a lower-field-traffic scalar control; all completeoperator results include reconciliation. S-vector uses first-touch, thread-private subdomains with duplicated boundaries and local reconciliation. E-vector stores every DoF per element for contiguous, conflict-free kernels, then reconciles copies to canonical DoFs. Both trade in-kernel indirection and conflicts for duplication and explicit synchronization. (A) Field locality. Starting from explicit-SIMD L-vector coloring, elastic SME at 𝑝 = 15 improves from 1663 GFLOP/s in L-vector to 1748 GFLOP/s in S-vector and 3447 GFLOP/s in kernel-only E-vector. Field reorganization raises absolute throughput for both backends rather than uniquely favoring SME. (B) Reconciliation cost. Including reconciliation reduces Evector to 1954 GFLOP/s: its 1.97× kernel advantage over S-vector narrows to 1.12× for the complete operator. Element-local storage therefore moves, rather than eliminates, much of the shared-DoF work.

primary coupled-vector target. Acoustic follows the same attenuation sequence, while its simpler pointwise map provides the scalar control; we quantify this contrast below. Compiler-vectorized. In Figure 7, elastic irregular reaches only 187–401 GFLOP/s with SVE and 182–470 GFLOP/s with SME, giving 0.97–1.17× SME/SVE—far below the 2.2× contraction proxy. Coloring is 1.09–3.59× faster than flat atomics and is therefore the L-vector baseline. Explicit SIMD. For acoustic, explicit-SIMD/no-SIMD performance spans 0.98–1.18× in regular geometry and 0.97–1.03× in irregular geometry, peaking at 1.18×. Its simpler pointwise map therefore shows substantially less mitigation than elastic. For elastic irregular, explicit SIMD improves coloring by 1.77–2.52× for SVE and 1.81–2.86× for SME (Figure 7). The high-order SME/SVE ratio rises from about 1.1× to 1.27×, but remains well below the contraction proxy. This distinction between absolute throughput and relative backend speedup is important: explicit SIMD removes a pointwise-physics bottleneck shared by both backends, but does not restore contraction to the dominant share of operator time. Residual attenuation. For elastic irregular, estimated contraction time is only 0.18–0.55× null time, while observed SME time is 1.07–1.22× null time (Figure 7). Irregular SoA also raises elastic null time by 1.4–2.0×; variable-field bandwidth is 80–129 GB/s and incremental coefficient bandwidth is 156–353 GB/s, both below the 400 GB/s STREAM reference (Figure 7). These values do not indicate saturated bulk bandwidth. Rather, coloring, indirect shared-DoF access, and fragmented coefficient streams prevent the operator from using STREAM-like delivery. Together, the overlap and bandwidth diagnostics identify structural field and coefficient movement—not contraction speed alone—as the remaining loss.

6.2

Data-Layout Recovery Across the Full Operator

After explicit SIMD exposes movement costs, we recover full-operator performance in two stages: first field locality, then coefficient delivery. Figure 9 summarizes this complete recovery sequence. Unless marked kernel-only, reported results include reconciliation.

6.1

Coefficient-Layout Recovery: AoSoV

After field-layout recovery, irregular geometry leaves coefficient delivery as the next attenuation layer. We evaluate this stage for both physics: panels C–D of Figure 9 present elastic results, while acoustic provides the scalar control. SoA stores each coefficient component as a separate elementmajor stream. AoSoV instead follows the vector-blocked AoSoA principle [14], grouping one SIMD width of GLL points with blockcontiguous components (Figure 10). Mixed AoSoV packs selected coefficients in FP16, then converts them for the same FP32 formulas and contraction backend. (C) Exposed coefficient penalty. With field access already element-local, switching from regular data to irregular SoA coefficients produces a 1.78× reconciliation-inclusive slowdown at 𝑝 = 15. This isolates metric and material streams as the next attenuation layer. Because the contiguous E-vector null path perturbs streaming, panel D uses the more stable S-vector null backend for bandwidth diagnosis. (D) Coefficient-traffic recovery. AoSoV raises effective coefficient delivery above the 400 GB/s STREAM reference for 𝑝 = 5–15; this overlap-aware metric is not literal sustained bandwidth. At 𝑝 = 15, reconciliation-inclusive E-vector performance relative to regular geometry improves from 0.56 with SoA to 0.78 with AoSoV and 0.89 with mixed AoSoV. The latter uses validated FP16 coefficient storage, not a general accuracy guarantee for arbitrary distorted meshes. For acoustic, the AoSoV/SoA effective-bandwidth ratio spans 0.91–1.35× across orders and reaches 1.35× at 𝑝 = 15, confirming that high-order vector-blocked coefficient delivery also benefits scalar stiffness. Thus field locality exposes the coefficient bottleneck, and vector-blocked coefficient storage recovers most of the remaining irregular-geometry penalty.

6.3

Best-Configuration Envelope and Roofline Placement

The elastic-irregular envelope spans L-vector, S-vector, and Evector fields; compiler-vectorized and explicit-SIMD pointwise physics; SoA, AoSoV, and mixed-AoSoV coefficients; and available execution and reconciliation plans. Invalid or tolerance-failing

Field-Layout Recovery: L-, S-, and E-Vectors

Figure 8 defines the L-vector/S-vector/E-vector semantics; panels A–B of Figure 9 present elastic results. Acoustic follows the same 7

Yinuo Wang, Lin Gan, Tianqi Mao, Zeyu Song, Wubing Wan, Jiayu Fu, Zekun Yin, Yuyang Jin, Xiaohui Duan, Wei Xue, and Guangwen Yang

28 31

29 32

30 33

34 36

35 37

38 40

39 41

42 44

17 20

18 21

22 24

23 25

41 44

Subdomain 0

Colored Elements

42 45

47 49

48 50

Subdomain 1

10 11

33 36

Full Result

Partial Result

(a) L-Vector

(b) S-Vector

Figure 8: Field storage and accumulation. L-vector uses shared DoFs with colored execution; S-vector and E-vector duplicate interface DoFs by subdomain and element, respectively. Storage indices remain contiguous across each rank-local vector. Partial interface results require reconciliation, whereas locally complete results do not. A. Field-Layout Recovery Under Regular Geometry

D. Coefficient-Layout Recovery

B. Reconciliation-Inclusive Time Decomposition

C. Irregular-Coefficient Penalty After Field Recovery

Figure 9: Field- and coefficient-layout recovery for elastic stiffness. (A) Regular-geometry L-vector, S-vector, and E-vector results isolate field locality. (B) Kernel and reconciliation-time decomposition shows where shared-DoF work moves. (C) Reintroducing irregular SoA coefficients isolates the remaining traffic penalty. (D) S-vector bandwidth diagnostics and reconciliation-inclusive E-vector results show recovery from AoSoV and mixed-AoSoV storage. results are excluded, and S-vector/E-vector times include reconciliation. Each backend independently selects its fastest remaining configuration at every order, so Figure 11 is a best-to-best rather than fixed-layout comparison. The envelope exposes the corresponding order-dependent transition. At 𝑝 = 4, best SME reaches only 0.94× best SVE, while even the null backend reaches only 1.23×. The short 𝑛 = 5 contraction occupies 5/16 of the FP32 matrix width and accounts for 60.4% of modeled FLOPs, so pointwise physics, coefficient loads and stores, field gathers and scatters, and reconciliation remain prominent. As

order increases, fuller matrix utilization and a rising contraction share generally raise SME/SVE to 1.64× at 𝑝 = 15, increasing the original irregular L-vector result from about 1.1× to about 1.6×. Yet null/SME remains only 1.11–1.32× across the sweep and is 1.11× at 𝑝 = 15. The null backend removes contraction instructions while preserving the surrounding operator path, so this ratio is both a diagnostic of non-contraction cost and an empirical counterfactual ceiling from making contractions free. Thus contraction FLOP share describes the increasing opportunity for matrix acceleration, whereas the null envelope measures the remaining elapsed-time 8

Beyond Fast Contractions: Attenuation and Recovery of Matrix-Engine Speedups in High-Order Finite Elements Component 1

Component 2

Component 3

1 2 3 4

Best validated local operator by backend Analytical TFLOP/s

Block 0

Component 0

FP16 Storage

Block 1

5 6 7 8

FP16 Load

SIMD Width 4 FP32 8 FP16

5 6 7 8 5 6 7 8 FP16 Vector

Vectors

SVE Shuffle Intrinsics Convert to Two FP32 Vectors

5 6 7 8

5 6 7 8 FP32 Vectors

Speedup over best SVE

Figure 10: AoSoV coefficient semantics and mixed-precision conversion. Each block stores one SIMD width of points contiguously by component. Mixed AoSoV packs selected components in FP16, then shuffles and converts each load into two FP32 vectors for pointwise computation.

4 2 0

00 M4

3.75 TFLOP/s

B/s

5 10 15 Arithmetic intensity (FLOP/byte)

Acoustic regular Acoustic irregular

Elastic regular Elastic irregular

SME SVE

Figure 12: Characteristic roofline after layout recovery, with intensity fixed to analytical work and the original traffic model. The horizontal compute roof is the 3.75 TFLOP/s SVE peak. Filled solid and open dashed curves denote the best validated configurations for SME and SVE, respectively.

1 0

4 SVE

7 9 11 13 15 Polynomial order p SME

Null

Discretization Implications

The recovered SME configurations attain their highest analytical FLOP rates at high order, but tensor-product stiffness performs 𝑂 (𝑛 4 ) work for 𝑂 (𝑛 3 ) unknowns, where 𝑛 = 𝑝 + 1. We therefore connect arithmetic rate to DoF throughput and an iso-accuracy cost projection in Figure 13. (1) Operator throughput. Reconciliation-inclusive GDoF/s rises toward high order for both backends (Figure 13). Improved hardware utilization therefore offsets the increasing work per DoF, showing that the high-order GFLOP/s gain is not merely an arithmeticcount artifact. (2) Iso-accuracy frontier. We consider dispersion-dominated error for a homogeneous scalar 20 Hz Ricker wave problem with Newmark time stepping. For receiver 𝑠, we compare complete observed and reference waveforms over [0,𝑇rec ] using ∫ 𝑇rec [𝑓obs,𝑠 (𝑡) − 𝑓ref,𝑠 (𝑡)] 2 d𝑡 𝜀𝑠 = 0 . (7) ∫ 𝑇rec 2 d𝑡 𝑓 (𝑡) ref,𝑠 0

Best SVE = 1

Figure 11: Best validated elastic-irregular speedup envelope. Each backend selects its fastest pointwise implementation, field/coefficient layout, and execution plan; bars are normalized to the best SVE configuration at each order.

We require the normalized squared waveform 𝐿2 error 𝜀𝑠 ≤ 10−2 . Group 2 first varies Δ𝑡 = 5×10−5 –8×10−4 at 𝑝 = 7 and ℎ = 40; we select Δ𝑡 acc = 2×10−4 where temporal error is comparable to spatial error. To isolate the spatial frontier, element-size sweeps at 𝑝 = 4, 7, 15 use Δ𝑡 = 10−4 , with ℎ ∈ {17.5, 20, 22.5, 25}, {40, 45, 50, 60}, and {110, 120, 125, 130, 135}, respectively. Error tables sample propagation distance from 0.5 to 12.49 km at 10 m intervals. For each sweep, station-wise error versus distance is fit by least squares to 𝑐 1𝑥 + 𝑐 2𝑥 2 . At each 𝑝, the frontier uses the largest tested ℎ whose fitted error at the receiver corresponding to approximately 10 s of propagation satisfies the criterion (Figure 13). This selects (𝑝, ℎ) = (4, 20), (7, 50), (15, 125). We use this homogeneous scalarwave calibration as a relative points-per-wavelength proxy for the elastic projection, assuming its order-to-order spatial-resolution trend transfers to smooth elastic propagation. Relative to 𝑝 = 4, the

benefit. This is not a strict hardware bound because removing contraction can perturb caching and overlap. Figure 12 places configurations selected under the same bestvalidated rule on the characteristic roofline [39], keeping intensity tied to analytical work and the original L-vector traffic model so layout recovery appears as vertical progress for both SVE and SME. Across acoustic and elastic, regular and irregular operators, recovered SME curves sit above the corresponding SVE curves at comparable intensity. Acoustic irregular SME reaches 80.5% of the SVE roof with FP32 AoSoV and 99.8% with mixed AoSoV at 𝑝 = 15, so its near-roofline result additionally depends on reduced coefficient precision. Elastic irregular SME reaches 56.2% and 64.1% of the SVE roof with FP32 and mixed AoSoV, respectively. Its qualitative recovery therefore does not depend on mixed precision, while remaining farther from the roof than the acoustic cases. 9

Yinuo Wang, Lin Gan, Tianqi Mao, Zeyu Song, Wubing Wan, Jiayu Fu, Zekun Yin, Yuyang Jin, Xiaohui Duan, Wei Xue, and Guangwen Yang

Group 1. Measured Operator Performance

Group 3. Iso-Accuracy Cost Projection

b) Throughput (GDoF/s)

a) Global DoF

b) Stable Timestep

c) Rel. Operator Time

a) Analytical GFLOPS

Measured Operator Cost

Accuracy Constraint

Group 2. Spatial–Temporal Error Calibration

Figure 13: From measured operator throughput to iso-accuracy stiffness cost. Group 1 reports reconciliation-inclusive elasticirregular GFLOP/s and GDoF/s. Group 2 uses station-wise waveform errors and fitted trends to determine the largest element size satisfying the normalized squared waveform 𝐿2 -error threshold 10−2 at each polynomial order. Group 3 combines global DoF count, stable time step, and measured operator time to expose the projected stiffness-cost tradeoff. latter points reduce the three-dimensional DoF count to 0.343 and 0.216, respectively, and also reduce pressure on the 4 GB-per-NUMA HBM capacity. (3) Projected stiffness cost. We combine the frontier with measured elastic-irregular throughput as follows, where Δ𝑡 CFL denotes the Courant–Friedrichs–Lewy (CFL) stability limit: 𝑇phys 𝑁 DoF proj 𝑇stiff ∝ . (8) min(Δ𝑡 acc, Δ𝑡 CFL ) 𝑅DoF

under-resolved features, and nonlinear stability can further require filtering or regularization. These effects are established concerns in high-order methods [2, 3, 22, 36, 38] and may shift the hardwarefavored frontier.

Related Work

Matrix Engines for Scientific Operators. Matrix engines now support a widening set of scientific operators beyond dense GEMM, including sparse–sparse multiplication [40], sparse–dense and sampled matrix products [15, 31], stencil computations [10, 17, 20, 26, 37, 41], and broader parallel patterns [27]. High-order FEM is a natural target because sum factorization exposes tensor contractions. LIBXSMM accelerates small CPU products [18]; GPU work optimizes tensor-product kernels [32] and maps them to tensor cores [11]. Concurrent work by Tu et al. programs FP64 tensor cores for fused full-operator MFEM kernels at scale, including a tsunami digital twin [19, 35].

Here 𝑅DoF is reconciliation-inclusive throughput and the step count is set by the smaller of the accuracy- and CFL-admissible time steps. Because minimum GLL spacing scales as 𝑂 (ℎ/𝑛 2 ), increasing ℎ competes with the high-order CFL penalty. The measured CFL limit decreases from 1.15 × 10−3 at 𝑝 = 4 to 1.05 × 10−3 at 𝑝 = 7 and 6.3 × 10−4 at 𝑝 = 15. The accuracy-selected step, 2 × 10−4 , is below all three limits, so projected stiffness speedup equals per-step speedup. From 𝑝 = 4 to 𝑝 = 7 it is 3.71× for SVE and 5.58× for SME. From 𝑝 = 7 to 𝑝 = 15 it is 1.84× and 2.14×. Under a CFL-limited sensitivity these become Limits of Prior Evidence and Our Position. Prior evidence consists 3.38× and 5.10×, then 1.10× and 1.28×. Thus 𝑝 = 7 is favorable largely of microkernel studies or implementation-to-implementation in either regime, whereas the added benefit of 𝑝 = 15 depends comparisons between matrix-unit and separate CUDA-core, SIMD, on time-step selection (Figure 13). These are stiffness-contribution or library paths. Reported gains may therefore combine matrix projections, not end-to-end speedups. instructions with changes in data reuse, work mapping, and kernel This frontier is deliberately a smooth, homogeneous wave-propagation fusion, making it difficult to attribute speedup specifically to the macase study rather than a universal prescription for high order. Realtrix unit. Coverage, data movement, and numerical format further istic simulations add geometric error from curved elements, heterogovern application benefit [13, 27]. We instead compare matched geneous material interfaces, and mesh-quality constraints; shocks, SVE, SME, and null backends on the same LX2 cores, then factor 10

Beyond Fast Contractions: Attenuation and Recovery of Matrix-Engine Speedups in High-Order Finite Elements

pointwise SIMD, field representation, reconciliation, and coefficient layout across the complete operator. Our contribution is an operator-level methodology rather than a new framework. Mainstream systems including libCEED [8], MFEM [4, 5], deal.II [7], Nek5000/NekRS [16], and SPECFEM3D [30] expose relevant matrix-free abstractions. Their optimized paths, however, do not always make the contraction backend, pointwise implementation, field representation, coefficient layout, and reconciliation strategy independently selectable. Our methodology identifies independent control of these interfaces as a key frameworklevel requirement for isolating matrix-engine benefit and recovering acceleration across the complete operator.

Methods Library. Computers & Mathematics with Applications 81 (2021), 42–74. doi:10.1016/j.camwa.2020.06.009 [5] Julian Andrej, Nabil Atallah, Jan-Phillip Bäcker, Jean-Sylvain Camier, Dylan Copeland, Veselin Dobrev, Yohann Dudouit, Tobias Duswald, Brendan Keith, Dohyun Kim, Tzanio Kolev, Boyan Lazarov, Ketan Mittal, Will Pazner, Socratis Petrides, Syun’ichi Shiraiwa, Mark Stowell, and Vladimir Tomov. 2024. HighPerformance Finite Elements with MFEM. The International Journal of High Performance Computing Applications (2024). doi:10.1177/10943420241261981 [6] Arm Ltd. 2022. Arm Architecture Reference Manual Supplement: The Scalable Matrix Extension. Arm Ltd. https://developer.arm.com/documentation/ddi0616/ latest [7] Daniel Arndt, Wolfgang Bangerth, Denis Davydov, Timo Heister, Luca Heltai, Martin Kronbichler, Matthias Maier, Jean-Paul Pelteret, Bruno Turcksin, and David Wells. 2021. The deal.II Finite Element Library: Design, Features, and Insights. Computers & Mathematics with Applications 81 (2021), 407–422. doi:10. 1016/j.camwa.2020.02.022 [8] Jed Brown, Ahmad Abdelfattah, Valeria Barra, Natalie Beams, Jean-Sylvain Camier, Veselin Dobrev, Yohann Dudouit, Leila Ghaffari, Tzanio Kolev, David Medina, Will Pazner, Thilina Ratnayaka, Jeremy Thompson, and Stan Tomov. 2021. libCEED: Fast Algebra for High-Order Element-Based Discretizations. Journal of Open Source Software 6, 63 (2021), 2945. doi:10.21105/joss.02945 [9] Claudio G. Canuto, M. Yousuff Hussaini, Alfio M. Quarteroni, and Thomas A. Zang. 2007. Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics. Springer-Verlag, Berlin, Heidelberg. [10] Yuetao Chen, Kun Li, Yuhao Wang, Donglin Bai, Lei Wang, Lingxiao Ma, Liang Yuan, Yunquan Zhang, Ting Cao, and Mao Yang. 2024. ConvStencil: Transform Stencil Computation to Matrix Multiplication on Tensor Cores. In Proceedings of the 29th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming. Association for Computing Machinery, 333–347. doi:10.1145/ 3627535.3638476 [11] Cu Cui. 2024. Acceleration of Tensor-Product Operations with Tensor Cores. ACM Transactions on Parallel Computing 11, 4, Article 34 (2024). doi:10.1145/ 3695466 [12] Michel O. Deville, Paul F. Fischer, and Ernest H. Mund. 2002. High-Order Methods for Incompressible Fluid Flow. Cambridge University Press. doi:10.1017/ CBO9780511546792 [13] Jens Domke, Emil Vatai, Aleksandr Drozd, Peng Chen, Yosuke Oyama, Lingqi Zhang, Shweta Salaria, Daichi Mukunoki, Artur Podobas, Mohamed Wahib, and Satoshi Matsuoka. 2021. Matrix Engines for High Performance Computing: A Paragon of Performance or Grasping at Straws?. In Proceedings of the 35th IEEE International Parallel and Distributed Processing Symposium. 1056–1065. doi:10.1109/IPDPS49936.2021.00114 [14] Xiaohui Duan, Kai Xu, Yuandong Chan, Christian Hundt, Bertil Schmidt, Pavan Balaji, and Weiguo Liu. 2017. S-Aligner: Ultrascalable Read Mapping on Sunway Taihu Light. In 2017 IEEE International Conference on Cluster Computing. IEEE, 36–46. doi:10.1109/CLUSTER.2017.16 [15] Ruibo Fan, Wei Wang, and Xiaowen Chu. 2024. DTC-SpMM: Bridging the Gap in Accelerating General Sparse Matrix Multiplication with Tensor Cores. In Proceedings of the 29th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 3. Association for Computing Machinery, 253–267. doi:10.1145/3620666.3651378 [16] Paul Fischer, Stefan Kerkemeier, Misun Min, Yu-Hsiang Lan, Malachi Phillips, Thilina Rathnayake, Elia Merzari, Ananias Tomboulides, Ali Karakus, Noel Chalmers, and Tim Warburton. 2022. NekRS, a GPU-Accelerated Spectral Element Navier–Stokes Solver. Parallel Comput. 114 (2022), 102982. doi:10.1016/j. parco.2022.102982 [17] Qiqi Gu, Chenpeng Wu, Heng Shi, Jianguo Yao, and Haibing Guan. 2026. Do We Need Tensor Cores for Stencil Computations? arXiv:2603.00477 [cs.DC] [18] Alexander Heinecke, Greg Henry, Maxwell Hutchinson, and Hans Pabst. 2016. LIBXSMM: Accelerating Small Matrix Multiplications by Runtime Code Generation. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE Press, Article 84, 11 pages. [19] Stefan Henneking, Sreeram Venkat, Veselin Dobrev, John Camier, Tzanio Kolev, Milinda Fernando, Alice-Agnes Gabriel, and Omar Ghattas. 2025. 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. Association for Computing Machinery, 60–71. doi:10.1145/3712285.3771787 [20] Han Huang, Jiabin Xie, Guangnan Feng, Xianwei Zhang, Dan Huang, Zhiguang Chen, and Yutong Lu. 2025. HStencil: Matrix–Vector Stencil Computation with Interleaved Outer Product and MLA. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. Association for Computing Machinery, 1816–1829. doi:10.1145/3712285.3759769 [21] Mark T. Jones and Paul E. Plassmann. 1993. A Parallel Graph Coloring Heuristic. SIAM Journal on Scientific Computing 14, 3 (1993), 654–669. doi:10.1137/0914041 [22] Andreas Klöckner, Tim Warburton, and Jan S. Hesthaven. 2011. Viscous Shock Capturing in a Time-Explicit Discontinuous Galerkin Method. Mathematical Modelling of Natural Phenomena 6, 3 (2011), 57–83. doi:10.1051/mmnp/20116303

Conclusion

Useful matrix-engine speedup is a property of the whole operator path. On LX2, the 4× SME peak advantage becomes about 2.2× for matched contractions and about 1.1× in the conventional irregular L-vector operator. Coordinated pointwise, field-, and coefficientlayout recovery raises the high-order elastic-irregular speedup to about 1.6×. Four operator-level lessons emerge. First, explicit SIMD for pointwise physics is necessary to expose backend differences that compilerdependent vectorization can mask. Second, the null backend is only 1.11–1.32× faster than recovered SME, providing an empirical ceiling on further contraction-only optimization. Third, continuousGalerkin shared-DoF cost is relocated rather than eliminated: Lvector, S-vector, and E-vector trade in-kernel indirection and conflicts against storage duplication and reconciliation. Finally, vectorblocked AoSoV aligns coefficient delivery with SIMD consumption; its FP32 recovery shows that the qualitative benefit does not depend on mixed precision. The measured ratios are LX2-specific, but the diagnostic is not: whenever a matrix unit accelerates only the contraction subset of an operator, pointwise work, field movement, coefficient delivery, and reconciliation jointly determine the realizable speedup. Exposing matrix units therefore requires full-operator representation contracts rather than an isolated fast contraction. The smooth wave-propagation case further suggests that recovered throughput can favor higher polynomial order, while extension to realistic geometries, heterogeneous media, and nonlinear problems remains future work.

References [1] Ahmad Abdelfattah, Valeria Barra, Natalie Beams, Ryan Bleile, Jed Brown, JeanSylvain Camier, Robert Carson, Noel Chalmers, Veselin Dobrev, Yohann Dudouit, Paul Fischer, Ali Karakus, Stefan Kerkemeier, Tzanio Kolev, Yu-Hsiang Lan, Elia Merzari, Misun Min, Malachi Phillips, Thilina Rathnayake, Robert Rieben, Thomas Stitt, Ananias Tomboulides, Stanimire Tomov, Vladimir Tomov, Arturo Vargas, Tim Warburton, and Kenneth Weiss. 2021. GPU Algorithms for Efficient Exascale Discretizations. arXiv:2109.05072 [cs.DC] [2] Mark Ainsworth. 2004. Discrete Dispersion Relation for ℎ𝑝 -Version Finite Element Approximation at High Wave Number. SIAM J. Numer. Anal. 42, 2 (2004), 553–575. doi:10.1137/S0036142903423460 [3] Mark Ainsworth and Hafiz Abdul Wajid. 2009. Dispersive and Dissipative Behavior of the Spectral Element Method. SIAM J. Numer. Anal. 47, 5 (2009), 3910–3937. doi:10.1137/080724976 [4] Robert Anderson, Julian Andrej, Andrew Barker, Jamie Bramwell, Jean-Sylvain Camier, Jakub Cerveny, Veselin Dobrev, Yohann Dudouit, Aaron Fisher, Tzanio Kolev, Will Pazner, Mark Stowell, Vassileios Tomov, Ido Akkerman, Johann Dahm, David Medina, and Stefano Zampini. 2021. MFEM: A Modular Finite Element 11

Yinuo Wang, Lin Gan, Tianqi Mao, Zeyu Song, Wubing Wan, Jiayu Fu, Zekun Yin, Yuyang Jin, Xiaohui Duan, Wei Xue, and Guangwen Yang

[41] Yiwei Zhang, Kun Li, Liang Yuan, Jiawen Cheng, Yunquan Zhang, Ting Cao, and Mao Yang. 2024. LoRAStencil: Low-Rank Adaptation of Stencil Computation on Tensor Cores. In Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis. IEEE Press, Article 53, 17 pages. doi:10.1109/SC41406.2024.00059

[23] Tzanio Kolev, Paul Fischer, Misun Min, Jack Dongarra, Jed Brown, Veselin Dobrev, Tim Warburton, Stanimire Tomov, Mark S. Shephard, Ahmad Abdelfattah, Valeria Barra, Natalie Beams, Jean-Sylvain Camier, Noel Chalmers, Yohann Dudouit, Ali Karakus, Ian Karlin, Stefan Kerkemeier, Yuhsiang Lan, David Medina, Elia Merzari, Aleks Obabko, Will Pazner, Thilina Rathnayake, Cameron W. Smith, Lukas Spies, Kasia Swirydowicz, Jeremy Thompson, Ananias Tomboulides, and Vassileios Tomov. 2021. Efficient Exascale Discretizations: High-Order Finite Element Methods. The International Journal of High Performance Computing Applications (2021). doi:10.1177/10943420211020803 [24] Dimitri Komatitsch, Gordon Erlebacher, Dominik Göddeke, and David Michéa. 2010. High-Order Finite-Element Seismic Wave Propagation Modeling with MPI on a Large GPU Cluster. J. Comput. Phys. 229 (2010), 7692–7714. doi:10.1016/j. jcp.2010.06.024 [25] Dimitri Komatitsch and Jeroen Tromp. 1999. Introduction to the Spectral Element Method for Three-Dimensional Seismic Wave Propagation. Geophysical Journal International 139, 3 (1999), 806–822. doi:10.1046/j.1365-246x.1999.00967.x [26] Xiaoyan Liu, Yi Liu, Hailong Yang, Jianjin Liao, Mingzhen Li, Zhongzhi Luan, and Depei Qian. 2022. Toward Accelerated Stencil Computation by Adapting Tensor Core Unit on GPU. In Proceedings of the 36th ACM International Conference on Supercomputing. Association for Computing Machinery, 1–12. doi:10.1145/ 3524059.3532392 [27] Yuechen Lu, Hongwei Zeng, Marc Casas, and Weifeng Liu. 2026. Characterizing Matrix Multiplication Units across General Parallel Patterns in Scientific Computing. In Proceedings of the 31st ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming. Association for Computing Machinery. doi:10.1145/3774934.3786456 [28] John D. McCalpin. 1995. Memory Bandwidth and Machine Balance in Current High Performance Computers. IEEE Technical Committee on Computer Architecture Newsletter (1995), 19–25. [29] Anthony T. Patera. 1984. A Spectral Element Method for Fluid Dynamics: Laminar Flow in a Channel Expansion. J. Comput. Phys. 54, 3 (1984), 468–488. doi:10. 1016/0021-9991(84)90128-1 [30] Daniel Peter, Dimitri Komatitsch, Yang Luo, Roland Martin, Nicolas Le Goff, Emanuele Casarotti, Pieyre Le Loher, Federica Magnoni, Qinya Liu, Céline Blitz, Tarje Nissen-Meyer, Piero Basini, and Jeroen Tromp. 2011. Forward and Adjoint Simulations of Seismic Wave Propagation on Fully Unstructured Hexahedral Meshes. Geophysical Journal International 186, 2 (2011), 721–739. doi:10.1111/j. 1365-246X.2011.05044.x [31] Jinliang Shi, Shigang Li, Youxuan Xu, Rongtian Fu, Xueying Wang, and Tong Wu. 2025. FlashSparse: Minimizing Computation Redundancy for Fast Sparse Matrix Multiplications on Tensor Cores. In Proceedings of the 30th ACM SIGPLAN Annual Symposium on Principles and Practice of Parallel Programming. Association for Computing Machinery, 312–325. doi:10.1145/3710848.3710858 [32] Kasia Swirydowicz, Noel Chalmers, Ali Karakus, and Tim Warburton. 2019. Acceleration of Tensor-Product Operations for High-Order Finite Element Methods. The International Journal of High Performance Computing Applications 33, 4 (2019), 735–757. doi:10.1177/1094342018816368 [33] TOP500. 2026. TOP500 List, June 2026. https://www.top500.org/lists/top500/ 2026/06/. Accessed July 2026. [34] Sara Touhami, Filippo Gatti, Fernando Lopez-Caballero, David Alejandro Castro Cruz, and Didier Clouteau. 2022. SEM3D: A 3D High-Fidelity Numerical Earthquake Simulator for Broadband (0–10 Hz) Seismic Response Prediction at a Regional Scale. Geosciences 12, 3 (2022), 112. doi:10.3390/geosciences12030112 [35] Jiqun Tu, Ian Karlin, John Camier, Veselin Dobrev, Tzanio Kolev, Stefan Henneking, and Omar Ghattas. 2026. Accelerating High-Order Finite Element Simulations at Extreme Scale with FP64 Tensor Cores. arXiv:2603.09038 [cs.DC] [36] Peter E. J. Vos, Spencer J. Sherwin, and Robert M. Kirby. 2010. From ℎ to 𝑝 Efficiently: Implementing Finite and Spectral/ℎ𝑝 Element Methods to Achieve Optimal Performance for Low- and High-Order Discretisations. J. Comput. Phys. 229, 13 (2010), 5161–5181. doi:10.1016/j.jcp.2010.03.031 [37] Yinuo Wang, Tianqi Mao, Lin Gan, Wubing Wan, Zeyu Song, Jiayu Fu, Lanke He, Wenqiang Wang, Zekun Yin, Wei Xue, and Guangwen Yang. 2026. SMEStencil: Optimizing High-Order Stencils on Arm Multicore Using SME Unit. IEEE Transactions on Parallel and Distributed Systems 37, 3 (2026), 651–665. doi:10.1109/TPDS.2025.3650515 [38] Z. J. Wang, Krzysztof Fidkowski, Rémi Abgrall, Francesco Bassi, Doru Caraeni, Andrew Cary, Herman Deconinck, Ralf Hartmann, Koen Hillewaert, H. T. Huynh, Norbert Kroll, Georg May, Per-Olof Persson, Bram van Leer, and Miguel Visbal. 2013. High-Order CFD Methods: Current Status and Perspective. International Journal for Numerical Methods in Fluids 72, 8 (2013), 811–845. doi:10.1002/fld.3767 [39] Samuel Williams, Andrew Waterman, and David Patterson. 2009. Roofline: An Insightful Visual Performance Model for Multicore Architectures. Commun. ACM 52, 4 (2009), 65–76. doi:10.1145/1498765.1498785 [40] Orestis Zachariadis, Nitin Satpute, Juan Gómez-Luna, and Joaquín Olivares. 2020. Accelerating Sparse Matrix–Matrix Multiplication with GPU Tensor Cores. Computers & Electrical Engineering 88 (2020), 106848. doi:10.1016/j.compeleceng. 2020.106848 12

Beyond Fast Contractions: Attenuation and Recovery of Matrix-Engine Speedups in High-Order Finite Elements

Analyze & explore related research in Conceptio

Related documents