---
title: 'MMStencil: Multicore & Matrix Stencil Optimization'
url: https://www.emergentmind.com/topics/mmstencil
type: topic
---

# MMStencil: Multicore & Matrix Stencil Optimization

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 [1310.8232]. 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 [2603.00477] [2310.16298] [2506.22969] [2507.11067]. 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" [1310.8232] 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$ [1310.8232].

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 [1310.8232].

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 [2603.00477]. "Stencil Matrixization" presents a matrixization strategy closely related to MMStencil, though built on vector outer products rather than dense matrix–matrix multiply [2310.16298]. "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 [2506.22969]. 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 [2507.11067]. 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 [1310.8232]. 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 [1310.8232].

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 [1310.8232]. 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$ [1310.8232].

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 $B_i$ such that $M_{\text{tile}} \le C_{\text{eff}}$, where $C_{\text{eff}}$ is the effective cache capacity available to the process [1310.8232]. On multicore sockets with shared caches, the paper recommended approximating $C_{\text{eff}}$ by dividing shared cache capacity by the number of cores sharing it, multiplied by a headroom factor $\gamma<1$ [1310.8232].

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 $L \approx 6r+1$ and stores were $S=1$, giving bytes $\approx s_e(L+S)$ for write-allocate stores when no non-temporal store was used [1310.8232]. For MPI domain decomposition in 3D with halo width $r$, the communicated bytes per time step per process were approximated by
$$
m \approx 2r \cdot s_e \cdot A \cdot (B_yB_z + B_xB_z + B_xB_y),
$$
with a latency–bandwidth communication model
$$
T_{\text{comm}} \approx \sum_{\text{faces}}(\alpha + \beta \, m_{\text{face}}).
$$
These relations positioned MMStencil within the standard tension between cache fit, halo surface-to-volume ratio, and multicore contention [1310.8232].

## 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 [1310.8232]. 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” [1310.8232].

MMStencil instead equalized cache usage across cores during selection. Its execution was divided into two phases. In the selection phase, a set of $nC$ block-size combinations was tested, each for $nI$ iterations, while candidate evaluation was coordinated across cores on the node [1310.8232]. In the verification phase, the selected block size was re-run for $nI$ iterations on all processes so that ActualTime could be compared with MinTime [1310.8232].

The three-stage coordinated selection proceeded as follows [1310.8232]:

1. **Stage 1 (per core)**: each process $P_{k,i}$ 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 $b^{\min}_{k,i}$.
2. **Stage 2 (per node)**: the per-core winners were aggregated using one of three target functions:
   - **Min-Min**
   - **Min-Worst (MWMB)**
   - **Min-Average**
3. **Stage 3 (global)**: the global best block $b_{\text{best}}$ was chosen by minimizing across nodes [1310.8232].

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 [1310.8232].

Candidate generation was defined by partitions per dimension. For a 3D domain of size $(S_x,S_y,S_z)$ with $P$ partitions per dimension, candidate block sizes per dimension were taken from
$$
b_l \in \left\{1,\frac{S_l}{P-1}+1,\frac{2S_l}{P-1}+1,\ldots,\frac{S_l}{P-1}(P-1)+1\right\},
$$
for $l \in \{x,y,z\}$, giving $nC=P^3$ candidates when the same $P$ was used in all dimensions [1310.8232]. Selection-phase cost scaled linearly in $nC$ and $nI$, and the paper explicitly described the tradeoff between larger search spaces and modest final runtime gains [1310.8232].

## 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 [1310.8232]. Three RTM problem sizes were studied: $200 \times 200 \times 800$, $300 \times 300 \times 800$, and $250 \times 250 \times 1500$ [1310.8232].

The algorithms compared were OB, MMMB, MWMB, and MAMB. Over 10 runs per case, the following averages were reported [1310.8232]:

| Problem size | OB average | MMMB average | MWMB average | MAMB average |
|---|---:|---:|---:|---:|
| $200 \times 200 \times 800$ | 963.2 s | 906.4 s | 631.8 s | 688.5 s |
| $300 \times 300 \times 800$ | 1648.2 s | 1537.0 s | 1295.0 s | 1443.8 s |
| $250 \times 250 \times 1500$ | 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 [1310.8232]. 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 [1310.8232]. 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 [1310.8232].

A central robustness claim was that MWMB consistently produced $\text{MinTime} \approx \text{ActualTime}$, with difference $<1\%$ [1310.8232]. The paper gave the following examples [1310.8232]:

| Problem size | MWMB block | MinTime | ActualTime | Total |
|---|---|---:|---:|---:|
| $200 \times 200 \times 800$ | $(15,15,143)$ | 4.00 s | 3.99 s | 628.8 s |
| $300 \times 300 \times 800$ | $(15,15,159)$ | 9.40 s | 9.35 s | 1295.1 s |
| $250 \times 250 \times 1500$ | $(47,15,495)$ | 13.60 s | 13.57 s | 1949.0 s |

By contrast, OB often exhibited large MinTime–ActualTime gaps, including a reported $+126.8\%$ divergence in the large case, which translated into much worse total times [1310.8232].

The search-space tradeoff was also quantified. In the medium case, increasing $nC$ 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 [1310.8232]. The paper therefore recommended matching $nC$ to expected amortization over repeated runs [1310.8232].

## 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$, radius $r$, element size $s_e$, and number of arrays $A$; 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 $\gamma \approx 0.6$ [1310.8232].

Candidate block sizes were to be constrained by the cache-fit relation
$$
M_{\text{tile}} = s_e \cdot A \cdot \prod_{i=1}^{d}(B_i+2r) \le C_{\text{eff}},
$$
with a preference for larger $B$ in the contiguous dimension and smaller $B$ in non-contiguous dimensions to limit random access and conflict [1310.8232]. The paper explicitly advised starting with small $P$, such as $P=5$ giving $nC=125$, and increasing to 10 or 20 only if the application would be rerun many times [1310.8232].

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 [1310.8232]. If sockets shared a last-level cache, per-socket tiles were to fit within the divided LLC capacity [1310.8232]. 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 [1310.8232].

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 [2603.00477]. In that setting, stencil updates are treated as localized contractions, while Tensor Cores implement fixed-shape matrix multiply-and-accumulate operations:
$$
D_{m \times n} = A_{m \times k} \times B_{k \times n} + C_{m \times n}.
$$
The paper emphasized two gaps between stencil and MMA: contraction mismatch and operand size requirements [2603.00477]. It identified flattening and decomposing as the two principal transformation families for making stencils “look like” MMA, and analyzed redundancy through a sparsity factor $\mathbb{S}$ and a fusion redundancy factor $\alpha$ [2603.00477]. 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 [2603.00477].

"Stencil Matrixization" took a different route by expressing stencil updates as sums of vector outer products accumulated into a stationary matrix register [2310.16298]. That paper explicitly positioned its method as close to MMStencil but distinguished it by using vector outer products rather than dense MMA [2310.16298]. It reformulated stencils in scatter mode, introduced “coefficient lines,” and derived a tile-level rank-1 decomposition
$$
B_{r:r+n-1,\,r:r+n-1} = \sum_{j=0}^{2r}\sum_{i=0}^{n+2r-1} u_{i,j} v_{i,j}^{T},
$$
thereby showing how a stencil operator restricted to a tile can be represented as a sum of outer-product masks [2310.16298]. On a simulated ARM SME-like platform, it reported speedups up to $4.71\times$ over vectorized baselines for some 2D and 3D stencils [2310.16298].

"SparStencil" generalized the same matrixization impulse to Sparse Tensor Cores with 2:4 structured sparsity [2506.22969]. 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 [2506.22969]. The work introduced Adaptive Layout Morphing, Structured Sparsity Conversion, and Automatic Kernel Generation, and reported up to $7.1\times$ speedup and $3.1\times$ on average over a state-of-the-art framework across 79 stencil kernels [2506.22969].

The 2025 paper titled "MMStencil: Optimizing High-order Stencils on Multicore CPU using Matrix Unit" extended the terminology directly into multicore CPU matrix units [2507.11067]. 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 [2507.11067]. 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 $2.1\times$, while enabling RTM applications to yield $1.8\times$ speedup versus a highly optimized industrial Nvidia A100 GPGPU version [2507.11067].

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 [1310.8232]. 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 [2603.00477] [2310.16298] [2506.22969].

Third, a 2025 paper uses MMStencil in its title for multicore CPU matrix-unit optimization of high-order stencils [2507.11067]. 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 [2507.11067]. 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 [2404.02218]—but those works concern compiler infrastructure rather than MMStencil as such. Likewise, out-of-core tiling in OPS [1709.02125] and geometric stencil selection in meshless finite differences [2202.06426] 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 [1310.8232]. On later GPUs and matrix-capable processors, the problem shifted toward matrixizing stencil operators to fit MMA, outer-product, or sparse-matrix engines [2603.00477] [2310.16298] [2506.22969] [2507.11067]. 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.

Source: https://www.emergentmind.com/topics/mmstencil