MMStencil: Multicore & Matrix Stencil Optimization

Updated 6 July 2026

MMStencil is a research strategy that refines stencil computations through node-aware block-size selection and later matrixization techniques for enhanced performance.
It employs a coordinated, multicore MPI profiling method that minimizes cache contention and balances workload across shared caches and NUMA architectures.
The approach expands to harness matrix engines by reformulating stencil updates into matrix multiply-accumulate operations, achieving significant speedups on modern GPUs and CPUs.

MMStencil is a name associated with multiple strands of stencil-computation research. In its most concrete and historically documented form, it denotes a node-aware methodology for choosing cache-efficient block sizes for multicore MPI stencil codes, introduced for a high-order Reverse Time Migration (RTM) kernel on multicore clusters (Sena et al., 2013). In later literature, the term is also used more broadly for stencil “matrixization” or “stencil-to-MMA” mappings that reformulate stencil updates to exploit matrix engines such as Tensor Cores, vector outer-product units, or matrix units on multicore CPUs (GU et al., 28 Feb 2026, Zhao et al., 2023, Li et al., 28 Jun 2025, Wang et al., 15 Jul 2025). The resulting terminological overlap makes MMStencil both a specific autotuning algorithm in multicore MPI optimization and, in a broader editorial sense, a family of architecture-aware stencil acceleration strategies.

1. Origins and scope

The 2013 work "Improving Memory Hierarchy Utilisation for Stencil Computations on Multicore Machines" (Sena et al., 2013) introduced MMStencil as a method for choosing cache-efficient block sizes for multicore MPI stencil codes. Its immediate target was regular-grid stencil computation in three dimensions, specifically a production RTM kernel used in seismic imaging. That kernel was described as a finite-difference PDE solver of 10th order in space and 2nd order in time, corresponding spatially to a 3D star stencil with radius $r=5$ (Sena et al., 2013).

In that formulation, MMStencil addressed a practical difficulty of migrating legacy MPI stencil applications from monocore clusters to multicore machines, where memory hierarchy is no longer dedicated to a sole core. The paper framed the central problem as the identification of an efficient block size that respects shared-cache behavior under multicore execution (Sena et al., 2013).

Subsequent literature uses the label more loosely. "Do We Need Tensor Cores for Stencil Computations?" discusses an “MMStencil-style stencil-to-MM mapping” on NVIDIA Tensor Cores (GU et al., 28 Feb 2026). "Stencil Matrixization" presents a matrixization strategy closely related to MMStencil, though built on vector outer products rather than dense matrix–matrix multiply (Zhao et al., 2023). "SparStencil" explicitly situates its contribution relative to what it describes as turning stencil updates into matrix multiply(-accumulate) so they can run on Tensor Cores or other matrix engines (Li et al., 28 Jun 2025). A 2025 paper titled "MMStencil: Optimizing High-order Stencils on Multicore CPU using Matrix Unit" states in its abstract that it introduces SIMD- and matrix-unit-based optimizations for 3D high-order stencils on multicore CPUs (Wang et al., 15 Jul 2025). This suggests that MMStencil has become polysemous: one meaning is a specific multicore MPI block-selection method, while another is a broader matrix-centric optimization idiom for stencil kernels.

2. Stencil model and computational setting

The original MMStencil work focused on regular-grid stencil computations in 3D and assumed standard row-major 3D arrays with contiguous memory in one dimension (Sena et al., 2013). Its RTM kernel had a working set much larger than cache and exhibited strong anisotropy: only one dimension, typically the innermost, was laid out contiguously in memory, so accesses in the other dimensions were expensive without tiling (Sena et al., 2013).

The paper discussed loop splitting and cache-aware tiling as principal optimizations. Loop splitting performed separate sweeps in $x$ , $y$ , and $z$ to reduce register pressure and enable SIMD, while cache-aware tiling selected 3D block sizes that better exploited the available memory hierarchy (Sena et al., 2013). Although the RTM kernel dominated the study, the methodology was said to generalize to common stencil patterns including the 2D 5-point stencil, the 3D 7-point stencil, the 3D 27-point stencil, and high-order 3D star stencils with $r>1$ (Sena et al., 2013).

The paper also gave a cache-fit model for a $d$ -dimensional star stencil with radius $r$ , element size $s_e$ , and $A$ arrays:

$M_{\text{tile}} = s_e \cdot A \cdot \prod_{i=1}^{d}(B_i+2r).$

The intended use was to choose tile extents $x$ 0 such that $x$ 1, where $x$ 2 is the effective cache capacity available to the process (Sena et al., 2013). On multicore sockets with shared caches, the paper recommended approximating $x$ 3 by dividing shared cache capacity by the number of cores sharing it, multiplied by a headroom factor $x$ 4 (Sena et al., 2013).

The same work described stencil arithmetic intensity and communication in halo-exchange form. For a single-array star stencil update, the approximate loads per update were $x$ 5 and stores were $x$ 6, giving bytes $x$ 7 for write-allocate stores when no non-temporal store was used (Sena et al., 2013). For MPI domain decomposition in 3D with halo width $x$ 8, the communicated bytes per time step per process were approximated by

$x$ 9

with a latency–bandwidth communication model

$y$ 0

These relations positioned MMStencil within the standard tension between cache fit, halo surface-to-volume ratio, and multicore contention (Sena et al., 2013).

3. Node-aware block selection in multicore MPI

The distinctive contribution of MMStencil in the 2013 sense was a coordinated, three-stage, node-aware profiling algorithm for block-size selection (Sena et al., 2013). The paper contrasted this method with naive dynamic profiling, called Original Block (OB), in which each MPI process tested different block sizes asynchronously on the same node. Because cores shared caches, this meant that different processes evaluated candidates under different cache contention levels. As a result, the block that appeared best during selection could underperform during subsequent verification, producing large divergence between “MinTime” and “ActualTime” (Sena et al., 2013).

MMStencil instead equalized cache usage across cores during selection. Its execution was divided into two phases. In the selection phase, a set of $y$ 1 block-size combinations was tested, each for $y$ 2 iterations, while candidate evaluation was coordinated across cores on the node (Sena et al., 2013). In the verification phase, the selected block size was re-run for $y$ 3 iterations on all processes so that ActualTime could be compared with MinTime (Sena et al., 2013).

The three-stage coordinated selection proceeded as follows (Sena et al., 2013):

Stage 1 (per core): each process $y$ 4 tested its assigned list of block candidates, with the same list and order across cores on the node, and selected its per-core best block $y$ 5.
Stage 2 (per node): the per-core winners were aggregated using one of three target functions:
- Min-Min
- Min-Worst (MWMB)
- Min-Average
Stage 3 (global): the global best block $y$ 6 was chosen by minimizing across nodes (Sena et al., 2013).

The Min-Worst variant, denoted MWMB, selected at Stage 2 the per-node block that performed well even in the worst core’s contention scenario. The paper reported that MWMB consistently delivered the best results across problem sizes and runs because it was robust to shared-cache pressure (Sena et al., 2013).

Candidate generation was defined by partitions per dimension. For a 3D domain of size $y$ 7 with $y$ 8 partitions per dimension, candidate block sizes per dimension were taken from

$y$ 9

for $z$ 0, giving $z$ 1 candidates when the same $z$ 2 was used in all dimensions (Sena et al., 2013). Selection-phase cost scaled linearly in $z$ 3 and $z$ 4, and the paper explicitly described the tradeoff between larger search spaces and modest final runtime gains (Sena et al., 2013).

4. Experimental evidence for the original MMStencil method

The 2013 experiments were conducted on the Oscar cluster with 40 BULL nodes, Gigabit Ethernet, and nodes containing two Intel Xeon E5430 2.66 GHz quad-core processors, 12 MB L2 cache per processor, and 16 GB RAM; experiments used up to 128 cores (Sena et al., 2013). Three RTM problem sizes were studied: $z$ 5, $z$ 6, and $z$ 7 (Sena et al., 2013).

The algorithms compared were OB, MMMB, MWMB, and MAMB. Over 10 runs per case, the following averages were reported (Sena et al., 2013):

Problem size	OB average	MMMB average	MWMB average	MAMB average
$z$ 8	963.2 s	906.4 s	631.8 s	688.5 s
$z$ 9	1648.2 s	1537.0 s	1295.0 s	1443.8 s
$r>1$ 0	3281.0 s	2439.7 s	1926.1 s	2043.7 s

For the small case, MWMB had very low variance with standard deviation 10.2 s, and its best run was 621.0 s versus 914.3 s for the best OB run (Sena et al., 2013). For the medium case, the best MWMB run was 1278.2 s versus 1279.4 s for the best OB run, but MWMB was more consistent on average and in worst-case behavior (Sena et al., 2013). For the large case, MWMB had the best average, though MAMB found one slightly faster block in one run; the paper nevertheless characterized MWMB as dominant overall (Sena et al., 2013).

A central robustness claim was that MWMB consistently produced $r>1$ 1, with difference $r>1$ 2 (Sena et al., 2013). The paper gave the following examples (Sena et al., 2013):

Problem size	MWMB block	MinTime	ActualTime	Total
$r>1$ 3	$r>1$ 4	4.00 s	3.99 s	628.8 s
$r>1$ 5	$r>1$ 6	9.40 s	9.35 s	1295.1 s
$r>1$ 7	$r>1$ 8	13.60 s	13.57 s	1949.0 s

By contrast, OB often exhibited large MinTime–ActualTime gaps, including a reported $r>1$ 9 divergence in the large case, which translated into much worse total times (Sena et al., 2013).

The search-space tradeoff was also quantified. In the medium case, increasing $d$ 0 from 125 to 1000 to 8000 raised selection time drastically—approximately 542 s, 3899 s, and 30428 s—while reducing final application runtime only from approximately 1375 s to 1294 s to 1203 s (Sena et al., 2013). The paper therefore recommended matching $d$ 1 to expected amortization over repeated runs (Sena et al., 2013).

5. Methodological principles and practitioner guidance

The 2013 paper articulated MMStencil as a practical workflow rather than solely an autotuning heuristic. Its recommendations included identifying stencil dimensionality $d$ 2, radius $d$ 3, element size $d$ 4, and number of arrays $d$ 5; choosing loop order so that the contiguous dimension is innermost; and estimating effective cache capacity per core using a headroom factor, with an initial recommendation of $d$ 6 (Sena et al., 2013).

Candidate block sizes were to be constrained by the cache-fit relation

$d$ 7

with a preference for larger $d$ 8 in the contiguous dimension and smaller $d$ 9 in non-contiguous dimensions to limit random access and conflict (Sena et al., 2013). The paper explicitly advised starting with small $r$ 0, such as $r$ 1 giving $r$ 2, and increasing to 10 or 20 only if the application would be rerun many times (Sena et al., 2013).

The multicore mapping assumptions were equally specific. The paper recommended one MPI rank per core, pinning ranks to cores, distributing ranks evenly across sockets on dual-socket nodes, and using first-touch initialization to place pages on the local NUMA node (Sena et al., 2013). If sockets shared a last-level cache, per-socket tiles were to fit within the divided LLC capacity (Sena et al., 2013). When performance jitter persisted, nearby block sizes could be tested to avoid cache set conflicts, and huge pages were suggested as a possible means of reducing TLB misses (Sena et al., 2013).

The broader implication is that MMStencil, in its original sense, is not merely a search routine but a coordinated methodology that integrates cache modeling, MPI rank placement, NUMA placement, and verification of robustness under realistic shared-cache contention. This suggests that the term’s enduring relevance derives from its attention to the mismatch between isolated microbenchmark timing and multicore production behavior.

6. Relation to matrixization and matrix-engine stencil research

Later stencil research reframed MMStencil as a matrixization concept. "Do We Need Tensor Cores for Stencil Computations?" described “MMStencil-style stencil-to-MM mapping” on NVIDIA Tensor Cores (GU et al., 28 Feb 2026). In that setting, stencil updates are treated as localized contractions, while Tensor Cores implement fixed-shape matrix multiply-and-accumulate operations:

$r$ 3

The paper emphasized two gaps between stencil and MMA: contraction mismatch and operand size requirements (GU et al., 28 Feb 2026). It identified flattening and decomposing as the two principal transformation families for making stencils “look like” MMA, and analyzed redundancy through a sparsity factor $r$ 4 and a fusion redundancy factor $r$ 5 (GU et al., 28 Feb 2026). Its key conclusion was that Tensor Core acceleration has a narrow “sweet spot” determined by roofline ridge points and by whether transformed kernels move from CUDA-core compute-bound operation into Tensor-Core memory-bound operation (GU et al., 28 Feb 2026).

"Stencil Matrixization" took a different route by expressing stencil updates as sums of vector outer products accumulated into a stationary matrix register (Zhao et al., 2023). That paper explicitly positioned its method as close to MMStencil but distinguished it by using vector outer products rather than dense MMA (Zhao et al., 2023). It reformulated stencils in scatter mode, introduced “coefficient lines,” and derived a tile-level rank-1 decomposition

$r$ 6

thereby showing how a stencil operator restricted to a tile can be represented as a sum of outer-product masks (Zhao et al., 2023). On a simulated ARM SME-like platform, it reported speedups up to $r$ 7 over vectorized baselines for some 2D and 3D stencils (Zhao et al., 2023).

"SparStencil" generalized the same matrixization impulse to Sparse Tensor Cores with 2:4 structured sparsity (Li et al., 28 Jun 2025). It explicitly described dense-MM stencil systems as masking sparsity with dense operands and proposed instead to transform stencil structure so that nonzeros aligned with the hardware’s 2:4 constraint (Li et al., 28 Jun 2025). The work introduced Adaptive Layout Morphing, Structured Sparsity Conversion, and Automatic Kernel Generation, and reported up to $r$ 8 speedup and $r$ 9 on average over a state-of-the-art framework across 79 stencil kernels (Li et al., 28 Jun 2025).

The 2025 paper titled "MMStencil: Optimizing High-order Stencils on Multicore CPU using Matrix Unit" extended the terminology directly into multicore CPU matrix units (Wang et al., 15 Jul 2025). According to its abstract, it analyzes matrix-based acceleration strategies for 3D high-order stencils, introduces SIMD and matrix-unit optimizations for strided memory accesses, alignment conflicts, and redundant accesses, proposes memory optimizations to boost on-package memory efficiency, and presents a multi-thread parallelism paradigm to address data-sharing challenges caused by the absence of shared data caches (Wang et al., 15 Jul 2025). The same abstract states that the method includes DMA-based inter-NUMA communication and that it outperforms state-of-the-art libraries on Nvidia A100 GPGPU by up to $s_e$ 0, while enabling RTM applications to yield $s_e$ 1 speedup versus a highly optimized industrial Nvidia A100 GPGPU version (Wang et al., 15 Jul 2025).

Taken together, these works indicate that MMStencil has acquired a second, broader meaning: a strategy of representing stencil computation in matrix-friendly form so as to exploit emerging matrix engines, whether dense, sparse, or outer-product based. A plausible implication is that the original 2013 emphasis on contention-aware block selection and the later matrixization literature address different layers of the same problem: one manages the memory hierarchy under multicore execution, while the other redefines the computational primitive to match specialized hardware.

7. Ambiguities, misconceptions, and historical position

A recurrent misconception is that MMStencil denotes a single universally recognized algorithm. The published record instead shows at least three distinct usages.

First, MMStencil is a specific 2013 methodology for block-size selection in multicore MPI stencil codes under shared-cache contention (Sena et al., 2013). This is the most fully documented meaning.

Second, the term is used descriptively in later literature for stencil-to-matrix-multiply mappings on accelerators. In this broader sense, MMStencil is not one implementation but a class of transformations that flatten, decompose, or otherwise matrixize stencil operators to suit Tensor Cores or analogous units (GU et al., 28 Feb 2026, Zhao et al., 2023, Li et al., 28 Jun 2025).

Third, a 2025 paper uses MMStencil in its title for multicore CPU matrix-unit optimization of high-order stencils (Wang et al., 15 Jul 2025). The abstract establishes the direction of that work, but the supplied detailed text does not enumerate its microkernels, scheduling, or empirical tables beyond the abstract claims (Wang et al., 15 Jul 2025). This means that the 2025 use is historically important but less documented in the presently available detail than the 2013 usage.

Another misconception is that MMStencil is synonymous with any stencil DSL or distributed-memory compiler framework. The literature does show adjacent developments—for example, a shared MLIR-based compilation stack for distributed-memory stencil DSLs spanning Devito, PSyclone, and the Open Earth Compiler (Bisbas et al., 2024)—but those works concern compiler infrastructure rather than MMStencil as such. Likewise, out-of-core tiling in OPS (Reguly et al., 2017) and geometric stencil selection in meshless finite differences (Davydov et al., 2022) address different meanings of “stencil” and should not be conflated with MMStencil’s regular-grid performance-optimization lineage.

Historically, the importance of MMStencil lies in how it tracks changes in hardware. On early multicore clusters, the problem was selecting block sizes robustly under shared caches and NUMA effects (Sena et al., 2013). On later GPUs and matrix-capable processors, the problem shifted toward matrixizing stencil operators to fit MMA, outer-product, or sparse-matrix engines (GU et al., 28 Feb 2026, Zhao et al., 2023, Li et al., 28 Jun 2025, Wang et al., 15 Jul 2025). This suggests continuity of purpose rather than continuity of implementation: across its usages, MMStencil denotes attempts to reconcile stencil regularity with the actual execution substrate of contemporary processors.