PaperlayerView exact revision on arXiv

Source-linked AI summary

A flexible algorithm for calculating pair interactions on SIMD architectures

Szilárd Páll, Berk Hess

arXiv:1306.1737v1physics.comp-ph

TL;DR

Traditional particle-pair algorithms struggle with memory access, data shuffling, and SIMD utilization. This paper groups particles into tunable spatial clusters, achieving 1.5 to 3 faster performance on 8-wide SIMD than traditional methods.

  • Problem

    Traditional particle-pair schemes suffer from irregular neighbor lists, scattered memory access, data shuffling, and limited SIMD efficiency.

  • Method

    The algorithm computes interactions between fixed-size spatial particle clusters, tuning cluster dimensions to SIMD width while reusing loaded data.

  • Results

    1.5 to 3 faster on 8-wide SIMD than traditional methods for typical atomistic molecular simulations.

  • Takeaways & Limitations

    The tunable cluster scheme maps efficiently to existing and future CPU and GPU architectures, including wider SIMD units.

  • Takeaways & Limitations

    The paper’s motivating standard implementations are limited by irregular, non-contiguous neighbor lists and memory-operation overheads.

Abstract

from arXiv · show

Calculating interactions or correlations between pairs of particles is typically the most time-consuming task in particle simulation or correlation analysis. Straightforward implementations using a double loop over particle pairs have traditionally worked well, especially since compilers usually do a good job of unrolling the inner loop. In order to reach high performance on modern CPU and accelerator architectures, single-instruction multiple-data (SIMD) parallelization has become essential. Avoiding memory bottlenecks is also increasingly important and requires reducing the ratio of memory to arithmetic operations. Moreover, when pairs only interact within a certain cut-off distance, good SIMD utilization can only be achieved by reordering input and output data, which quickly becomes a limiting factor. Here we present an algorithm for SIMD parallelization based on grouping a fixed number of particles, e.g. 2, 4, or 8, into spatial clusters. Calculating all interactions between particles in a pair of such clusters improves data reuse compared to the traditional scheme and results in a more efficient SIMD parallelization. Adjusting the cluster size allows the algorithm to map to SIMD units of various widths. This flexibility not only enables fast and efficient implementation on current CPUs and accelerator architectures like GPUs or Intel MIC, but it also makes the algorithm future-proof. We present the algorithm with an application to molecular dynamics simulations, where we can also make use of the effective buffering the method introduces.

1. Introduction

The introduction motivates SIMD-friendly pair-interaction algorithms by identifying irregular, non-contiguous neighbor lists and costly memory access as limitations of traditional particle-wise schemes. It proposes fixed-size particle clusters to address these issues while accepting extra zero-valued interactions, and describes implementations across CPUs, GPUs, and molecular-dynamics software.

  • Motivation: Pair interactions with limited spatial range consume more than half of computational time in most particle simulations.Pair-correlation analysis also involves evaluating limited-range particle interactions.
  • Limitations of existing methods: Traditional double-loop neighbor-list algorithms reduce quadratic cost to linear by exploiting spherical cut-offs, but their neighbor lists are irregular and non-contiguous for data parallelization.Verlet lists and linked-cell algorithms are widely used, particularly in molecular-dynamics codes.
  • Architectural challenges: Cut-offs make neighbor data non-sequential, increasing shuffle costs on SIMD CPUs and causing scattered accesses, bandwidth waste, and memory-bound execution on GPUs.The relative shuffle cost depends on the cost of one pair interaction and the SIMD width.
  • Proposed method: The proposed approach computes interactions between fixed-size clusters of particles to regularize data parallelization and address irregular neighbor-list sizes and non-contiguous layouts.It may calculate extra interactions outside the cut-off, which evaluate to zero, but the performance gain outweighs this extra cost.
  • Implementation and scope: The algorithm operates at the interaction-kernel level, supports existing parallelization strategies, and was implemented for x86 SIMD architectures, NVIDIA GPUs, and GROMACS molecular-dynamics simulations.Implemented kernels cover SSE2, SSE4.1, AVX, AVX+FMA, and several Lennard-Jones and electrostatic interaction forms.

2. The algorithm

The algorithm groups particles into fixed-size spatial clusters that serve as computational units mapped to SIMD hardware. Cluster-based pair lists increase data reuse and reduce memory operations, while allowing tunable cluster sizes and requiring compact spatial grouping despite extra included pairs.

  • Cluster-based SIMD execution: The cluster algorithm increases data reuse by a factor of M, with the standard Verlet-list implementation represented by M=1 and N=1.The algorithm loads a cluster of M particles and loops over neighboring clusters of N particles.
  • Cluster-based SIMD execution: Fixed-size spatial clusters are mapped directly to SIMD units, with cluster sizes tuned to the hardware width.For a 4-way SIMD unit, particles are grouped in clusters of 4.
  • Cluster-based SIMD execution: A 4×4 cluster setup computes 16 pair interactions while loading and storing only 4 j-particles.Compared with the 1×1 scheme, this computes 4 times as many interactions per particle load/store and reduces memory operations.
  • Cluster pair-list construction: Pair-list data are represented as cluster pairs rather than individual particle neighbors, and Newton’s third law stores each pair only once.The cluster pair list is a generalized neighbor list; the classical neighbor list is the 1×1 case.
  • Cluster pair-list construction: Compact fixed-size clusters are generated by spatial gridding and binning, using grid spacing of (max(M, N)/ρ)1/3.Bounding boxes then support efficient cluster-pair construction through distance checks for M × N particle pairs.
  • Cluster pair-list construction: Cluster pair lists inherently include extra particle pairs, whose fraction rises rapidly with cluster size and falls rapidly with cut-off radius.The efficiency gains are intended to outweigh the cost of calculating these additional interactions.

3. Pair interaction kernels

The pair-interaction kernels are optimized by tuning cluster dimensions, masking, and SIMD implementation choices to keep computational units busy while limiting memory and latency costs. CPU and CUDA kernels use architecture-specific strategies, including M=4 on CPUs and M=8, N=4 on GPUs, with efficient Lennard-Jones and electrostatic parameter handling.

  • Implementation choices: CPU kernels use C with SSE and AVX intrinsics, while GPU kernels target NVIDIA CUDA because its tools are mature and performant.Modern compilers typically optimize intrinsic-based CPU code better across architectures than equivalent hand-written assembly.
  • Kernel design: Masking is usually more efficient than conditionals for exclusion and cut-off checks, using CPU bitwise AND operations and GPU multiplication by 0 or 1.GPU cut-off checks still use a conditional, which can reduce instruction efficiency.
  • CPU kernels: M=4 gives the best CPU performance because increasing M raises the compute-to-memory ratio by reusing each loaded j-particle across M interactions.The kernels aim to avoid stalls from memory-operation dependencies and instruction latencies.
  • GPU kernels: CUDA kernels are instruction-latency limited and use M=8, N=4 so 32 warp threads calculate an entire cluster-pair interaction simultaneously.Choosing N smaller than M increases computation per memory operation, while tight pair-list packing is required.
  • Interaction models: Geometric or Lorentz-Berthelot combination rules reduce costly SIMD parameter loading, while smooth PME corrections support linear table interpolation at full single precision.Arbitrary pair-parameter loads require many load and shuffle operations; PME interpolation can use a limited table size.

4. Performance in practice

The evaluation compares optimized kernel variants primarily on Intel Sandy Bridge, using cycle-based performance metrics and separate CPU/GPU compilation setups. Results show that wider SIMD and larger cluster configurations improve pair throughput, while shuffle and data-movement overheads limit gains.

  • Benchmark setup: Performance is evaluated on Intel Sandy Bridge, whose AVX support enables direct comparison of 4×4 versus 4×8 and 128-bit versus 256-bit setups.Performance data are reported in cycles, depending on microarchitecture rather than the exact CPU model.
  • Performance metrics: The evaluation reports pairs/cycle, IPC, flops/pair, and flops/cycle to measure raw performance, hardware utilization, arithmetic intensity, and computational throughput.Pairs/cycle is identified as the relevant raw-performance measure; flops/pair is minimized while flops/cycle is maximized.
  • CPU performance: The 256-bit RF kernel is only 13% faster than the 128-bit variant despite similar IPC, because shuffle and data-load overheads limit the benefit of wider SIMD.The kernels execute the same arithmetic instructions, so the small gain is attributed to data movement rather than computation.
  • Cluster-size comparison: The 4×2 RF kernel is 26% faster than the 2×2 kernel, while the 4×4 scheme delivers 50% higher performance than 4×2 in 256-bit kernels.The results rule out M = 2 as a viable option and support considering M = 4 or larger.
  • Architecture comparison: Sandy Bridge with 256-bit SIMD provides a 20% higher pair rate than 128-bit SIMD on Sandy Bridge, whereas Bulldozer is only marginally faster than 128-bit Sandy Bridge SIMD.The passage also reports significantly higher CUDA GPU performance when comparing one streaming multiprocessor with one CPU core.

5. Effective pair list buffering

Effective pair-list buffering evaluates performance using only non-zero interactions within the cut-off and can exploit extra cluster-scheme pairs in molecular dynamics. The 4×4 and 8×4 schemes support update frequencies around 10 steps, while optimized 1×1 performance favors roughly 10–15 steps.

  • Buffering mechanism: Molecular dynamics can exploit extra cluster-scheme pairs through Verlet-buffered pair lists, which permit particle motion without immediately invalidating the list.The buffer ensures interacting pairs are not missed while allowing some movement between list updates.
  • Buffering assumptions: Small pair-list imperfections can be tolerated in constant-temperature simulations because a thermostat removes excess heat from energy drift.The passage frames this as a practical alternative to requiring perfect pair-list conservation.
  • Performance measure: Effective performance counts only interactions within the cut-off radius calculated per cycle, excluding extra zero-valued buffer and cluster-pair interactions.This measure is intended to reflect absolute performance rather than the total number of calculated pairs.
  • Update frequency: The 1×1 pair-list construction takes four times longer than one interaction calculation, whereas 4×4 search and interaction calculation take about equal time.The 1×1 search is not fully optimized, and making it twice as fast would place its optimal update frequency between 10 and 15 steps.
  • Update frequency: The optimal pair-list update frequency is around 10 steps for both 4×4 and 8×4 schemes.For an optimized 1×1 search, the estimated optimum is between 10 and 15 steps.

6. Conclusions

The paper concludes that traditional SIMD pair-interaction algorithms are limited by memory, shuffling, scheduling, and pair-list constraints. The proposed method performs 1.5 to 3 times faster for typical atomistic simulations on 8-wide SIMD, exceeds 60% of peak flop rate on Sandy Bridge CPUs and CUDA GPUs, and maps to current and future architectures, though extra zero interactions can offset gains for expensive interactions.

  • Conclusions: Traditional particle-based SIMD algorithms are limited by high memory-to-arithmetic ratios, data shuffles, instruction scheduling, and pair-list construction constraints.These limitations reduce memory-latency hiding and have reached the limits of standard loop-unrolled approaches.
  • Conclusions: Extra zero interactions can outweigh the gains for expensive interactions.This is a stated limitation of the method’s setup.
  • Conclusions: 1.5 to 3 faster on 8-wide SIMD is achieved for typical atomistic molecular simulations compared with traditional methods.The stated comparison is for the proposed method versus traditional methods.
  • Conclusions: Above 60% of peak flop rate is achieved on Intel Sandy Bridge CPUs and CUDA GPUs.This performance is reported for the proposed method in typical atomistic molecular simulations.
  • Conclusions: The scheme maps well to future CPU and GPU architectures as well as existing architectures not discussed in the paper.Its floating-point operations per load/store operation can be tuned, enabling adjustment of arithmetic cycles per kernel.

Appendix

The appendix derives an upper bound on average energy drift from finite Verlet buffers using Gaussian displacement distributions, while noting that interactions narrow these distributions and make the bound conservative.

  • Energy-drift bound: The displacement variance for a freely moving particle is σ2 = t kBT/m, while the distance variance for two non-interacting particles is σ2 = t kBT (1/m1 + 1/m2).These Gaussian distributions have zero mean; particle interactions narrow the actual displacement distribution, so ignoring them gives an upper bound.
  • Energy-drift bound: The method uses a non-bonded cut-off rc and pair-list cut-off rℓ= rc + rb, where rb is the Verlet buffer size.The finite buffer size enters the energy-drift bound through the separation between the interaction and pair-list cut-offs.
  • Approximation: Small energy drift requires σ to remain small compared with both rc and rℓ, making the Gaussian-tail approximations accurate.The approximation is justified because the Gaussian distribution decays rapidly when σ is small relative to the cut-offs.
  • Total drift: Total energy drift is obtained by averaging over all particle pairs and weighting by particle count.The Gaussian function G(x) has zero mean and unit variance, with E(x) = 1.
Loading 1306.1737v1…