MatrixPIC: Accelerating PIC Simulations

• MatrixPIC is a high-performance framework that reformulates PIC current deposition as a matrix outer-product to enhance computational efficiency and memory throughput.
• The hybrid MPU–VPU execution model and GPMA-based incremental sorting optimize data vectorization and spatial locality, thereby reducing register traffic and latency.
• Performance evaluations show up to 8.7× kernel speedup and 83% of FP64 peak efficiency, demonstrating significant advancements over traditional PIC implementations.

MatrixPIC refers to both a cryptosystem for key exchange and, more prominently in current computational science, a high-performance algorithmic and architectural framework for accelerating Particle-in-Cell (PIC) simulations on emerging CPU hardware featuring Matrix Processing Units (MPUs). The following account centers strictly on the MatrixPIC system as precisely detailed in "Matrix-PIC: Harnessing Matrix Outer-product for High-Performance Particle-in-Cell Simulations" (Rao et al., 13 Jan 2026), a work introducing a matrix-centric redesign of the current deposition kernel in fully kinetic plasma simulation codes. MatrixPIC embodies a tightly co-designed approach, integrating block-matrix numerics, a hybrid SIMD/VPU–MPU execution pipeline, optimized memory layout, and an amortized O(1) incremental sorting strategy.

1. Mathematical Reformulation of PIC Deposition

MatrixPIC fundamentally recasts the standard "scatter-add" operation for current deposition in Particle-in-Cell frameworks. Traditional deposition equations, where each particle’s charge $q_p$ and velocity $\mathbf{v}_p$ updates the grid current at the eight nearest nodes using a shape function $S(x_i-x_p)$, are factorized (for the first-order 3D case) into coordinated 1D weight arrays $s_{x,i}(\Delta x_p)$, $s_{y,j}(\Delta y_p)$, $s_{z,k}(\Delta z_p)$ for $i,j,k\in\{0,1\}$.

The key mathematical transformation is grouping two particles and expressing their contributions to the eight-nodal "rhocell" buffer as a $4\times8$ matrix outer product:

• For particles $p_1$, $p_2$:
  • Construct vector $a\in\mathbb{R}^4$: stacked $s_{x,i} \cdot W_p^x$ for both particles.
  • Construct vector $b\in\mathbb{R}^8$: outer product of $s_{y,j}s_{z,k}$ terms, concatenated for both particles.
  • Their interaction yields a $C=a\otimes b\in\mathbb{R}^{4\times8}$ matrix whose sub-blocks correspond directly to the contributions of each particle.

Generalization to higher-order (e.g., third-order QSP) schemes increases block size, producing $m=n=8$ matrices ($64$ updates per pair). The total deposition over a particle tile is written as a high-density rank-1 update: $RhocellBlock += A \cdot B^\top$.

This block-matrix algebra maps precisely onto the hardware MPU’s outer-product-accumulate (MOPA) primitives, enabling bulk updates to memory-resident grid tiles.

2. Hybrid MPU–VPU Execution Model

MatrixPIC implements a synergistic pipeline integrating both classic wide SIMD (VPU) and the novel MOPA-capable MPUs:

• VPU phase: Data for particle tiles (positions, velocities, charges, weights) are loaded and vectorized. The VPU computes cell indices, intra-cell coordinates, 1D shape factor arrays, and forms particle move lists.
• Incremental sorting (see Section 3) maintains semi-sorted order within particle tiles.
• MPU phase: For each cell (in tile-sorted order), batches of two particles are processed. VPU registers assemble $A, B$ matrices; MPOA instructions issue $C_\text{tile}=A\otimes B$ and deposit results into the rhocell buffers.
• Final reduction: VPU performs a reduction, scalarly adding rhocell buffers into the global current grid.

By maximizing sequential cell access and batch update locality, the design drastically limits costly VPU–MPU register traffic, exposes high arithmetic intensity to the MPU, and minimizes memory synchronization across threads.

3. Amortized O(1) Incremental Sorting via GPMA

To maintain spatial locality (critical for MPU execution efficiency and cache utilization), MatrixPIC employs a Gapped Packed-Memory Array (GPMA) per particle tile:

• The GPMA structures particle indices in a contiguous array interspersed with preallocated gaps (“empty slots”).
• As particles move between cells due to motion (bounded by the CFL condition), insertions and deletions typically proceed in O(1) time by finding or exploiting an adjacent gap.
• Only when a cell block’s density crosses a threshold does the entire tile undergo a redistribution to re-establish even gap placement (worst-case O($N_\text{tile}$), but amortized to O(1) per operation).
• Algorithmically, this involves iterating through the pending move list, attempting to fill or vacate slots, shifting minimal data, and very rarely rebuilding the entire array.

This mechanism supports efficient, nearly contention-free locality maintenance without burdening the critical path for typical particle motions.

4. Memory and Data-Layout Co-Design

MatrixPIC couples its algorithmic advances to a rigorously tuned memory layout:

• Particles are stored as Structure-of-Arrays (SoA) blocks, tile-sized at $8\times8\times8$ to match the VPU width and maximize tile-level L1 cache residency.
• Each rhocell is a $[8\times N_{\text{cells}}]$ 2D array, where each column is 8 contiguous doubles (aligned with 64-byte cache lines), supporting direct MPU register-level writes and perfectly unit-stride VPU reductions.
• This layout forms the substrate for efficient VPU/MPU cooperativity—allowing vectorized, aligned loading, minimal cache thrashing, and maximal MOPA utilization.

Domain decomposition into particle tiles further ensures that all data required per update batch remains within the working set of relevant CPU caches.

5. Performance Evaluation on Emerging HPC Hardware

MatrixPIC was evaluated on a dual-die LX2 CPU platform (256 cores, each with 512-bit VPU and 8×8 FP64 MPU at 1.3 GHz) and contrasted with an NVIDIA A800 GPU:

• First-order CIC with 128 particles per cell (PPC): End-to-end runtime reduced from 74.13 s (WarpX baseline) to 24.90 s—a $2.98\times$ total speedup, and the compute phase (kernel) improved by $3.76\times$.
• Third-order QSP scheme: Kernel accelerated by $8.7\times$ versus baseline and $2.0\times$ over the fastest prior VPU implementation.
• Peak efficiency: MatrixPIC approaches $83.08\%$ of FP64 theoretical peak on LX2, which is $2.8\times$ that of a hand-tuned CUDA kernel on the NVIDIA A800 GPU.
• LWFA workloads: Up to $2.63\times$ total-time speedup, with benefits increasing in denser regimes ($\text{PPC}\geq 32$).

Performance tables confirm both wall-time and efficiency gains, with particular superiority in scenarios where MPU block-outer-products can be heavily exploited.

6. Design Constraints, Trade-Offs, and Extensibility

MatrixPIC’s architectural and algorithmic trade-offs include:

• At low particle densities ($\text{PPC}<8$), the sorting and VPU logic overheads are not amortized; the framework is designed to fallback to more suitable scalar/VPU implementations in such regions.
• The current method addresses direct (non-charge-conserving) deposition only. Adapting Esirkepov or Vay charge-conserving methods within a matrix-centric paradigm requires extended continuity-preserving kernel reforms.
• Hardware specificity: The framework fundamentally depends on availability of MOPA-capable MPUs (such as Arm SME, Apple M4, NVIDIA Grace, AMD MI300) and is not portable to legacy SIMD-only CPUs without loss of advantage.
• Occasional GPMA rebuilds trade small, infrequent latency for O(1) amortized performance, consistent with sorted-tile strategies in modern kinetic codes.

Potential extension vectors include gather-field interpolation via outer-product mapping, auto-tuning tile and block sizes for adaptive density, and migration to new MPU-rich architectures or, if scatter primitives are feasible, to tensor-core style GPU accelerators.

7. Impact and Implications for PIC and Computational Plasma Physics

MatrixPIC demonstrates that classical, sparsely updated PIC kernels can be dramatically accelerated by reformulating them as dense, block-outer-product operations and exploiting hardware innovation in matrix processing. By aligning algorithmic structure, particle sorting, and data layout to the strengths of modern CPU MPUs, MatrixPIC achieves kernel speedups in the $2.98\times$–$8.7\times$ range and over $80\%$ of peak floating-point throughput on leading-edge CPU platforms. This approach substantially outperforms both hand-optimized SIMD CPU and CUDA GPU rivals in mainstream deposition workloads (Rao et al., 13 Jan 2026).

This advancement suggests a paradigmatic shift for simulation software: by co-designing numerical kernels for hardware primitives (particularly outer-product and collective accumulation), substantial gains are realizable, especially in high-density, cache-sensitive, or chronically contention-limited applications. In computational plasma physics, MatrixPIC provides a template for future high-performance PIC engines targeting exascale and beyond.