PaperlayerView exact revision on arXiv

Source-linked AI summary

Performance Portable $\mathrm{SU}(N)$ Lattice Gauge Theory Simulation with Kokkos

Wei Sun

arXiv:2610.09766v1hep-latcs.DChep-phphysics.comp-ph

TL;DR

Separate architecture-specific implementations make portable lattice gauge simulation costly to maintain. kwqft implements Wilson pure-gauge Monte Carlo in Kokkos for compile-time SU(N) and D, and its verification and performance tests report competitive results across multiple platforms.

  • Problem

    Heterogeneous hardware makes separate architecture-specific lattice-gauge implementations costly to maintain, while portable source alone does not ensure efficient execution.

    Source: §1 · p. 1

  • Method

    kwqft implements Wilson pure-gauge Monte Carlo with compile-time SU(N) and dimension parameters, Kokkos kernels, host SIMD, and MPI halo exchange overlapped with interior updates.

    Source: §1 · p. 2, §7 · pp. 12–13

  • Results

    Across verification and portability tests, kwqft reproduces reference plaquettes and performs competitively across platforms; strong scaling of a 128^3 × 256 SU(4) lattice from 16 to 512 nodes reaches 16.2× speedup.

    Source: §4.2 · p. 7, §6.5 · p. 12, §7 · p. 13

  • Takeaways & Limitations

    The results support writing a complete pure-gauge generation and measurement path once in Kokkos for the evaluated CPU, GPU, and distributed-memory settings.

    Source: §6.5 · p. 12, §7 · p. 13

    The implementation is currently pure gauge; extending it to dynamical fermions remains a longer-term objective requiring additional fields, solvers, force kernels, integrators, and an HMC algorithm.

    Source: §7 · p. 13

Abstract

from arXiv · show

The increasing diversity of high performance computing systems makes separate, architecture specific implementations of lattice gauge theory algorithms costly to maintain. We present \texttt{kwqft}, a performance portable Kokkos implementation of Wilson pure gauge Monte Carlo simulation for $\mathrm{SU}(N)$ Yang-Mills theory in an arbitrary number of space-time dimensions. The gauge group order $N$ and the dimension $D$ are compile time parameters. A single source targets the Serial, OpenMP, CUDA, HIP, and SYCL execution spaces, with MPI halo exchange overlapped with interior updates. The implementation reproduces the exact two-dimensional plaquette and published three and four dimensional values for gauge groups up to $\mathrm{SU}(17)$. On an NVIDIA A100 the Kokkos CUDA backend is competitive with a native CUDA code, SIMD acceleration improves the OpenMP path on Armv9 processors, and a large scale speedup is demonstrated for an $\mathrm{SU}(4)$ lattice on the LineShine supercomputer, currently ranked first on the TOP500 list.

1. Introduction

Heterogeneous hardware makes separate architecture-specific implementations costly to maintain, while portable source alone does not ensure efficient execution. The paper presents kwqft as a Kokkos pure-gauge implementation designed for general SU(N), dimension-generic simulation and competitive performance across platforms.

1. Introduction

  • Heterogeneous CPUs and accelerators make maintaining separate CUDA, OpenMP, HIP, and SYCL implementations costly, and compilation alone does not guarantee efficient execution.Performance portability seeks one implementation that serves declared platforms without sacrificing an unacceptable fraction of their performance.Source: §1 · p. 1
  • Earlier portable lattice-QCD studies focused primarily on fermion operators or SU(3), leaving the complete pure-gauge path for general N and D as the question addressed here.The paper also contrasts its goals with CUDA-specific, fixed-four-dimensional pure-gauge work.Source: §1 · pp. 1–2
  • kwqft combines compile-time SU(N) and dimension-generic updates, Kokkos execution, host SIMD, and MPI halo exchange overlapped with interior updates.The contribution list also identifies physics verification and performance studies across an A100, Armv9, and up to 512 LineShine nodes.Source: §1 · p. 2

2. Lattice formulation and update algorithm

The formulation uses the Wilson gauge action and a local staple to update links. Pseudo-heatbath, overrelaxation, checkerboard parallelism, and standard observables complete the simulation trajectory.

2.1. Wilson pure gauge action

  • The paper focuses on the isotropic Wilson action, with anisotropic spatial and temporal plaquettes assigned weights β/ξ0 and βξ0, respectively.The anisotropy parameter is ξ0, and the isotropic case sets ξ0 = 1.Source: §2.1 · p. 2
  • For a selected link, the staple collects the action-dependent neighboring-link products and is the update algorithms’ only stencil operation.It involves 2(D−1) three-link products and reads 6(D−1) distinct links.Source: §2.1 · p. 3

2.2. Cabibbo–Marinari pseudo-heatbath

  • For N > 2, Cabibbo–Marinari pseudo-heatbath applies SU(2) updates across all color-pair subgroups; SU(2) instead has a direct sampler.The SU(2) draw uses Kennedy–Pendleton at large k and Creutz at small k.Source: §2.2 · p. 3

2.3. Overrelaxation

  • Microcanonical overrelaxation reflects each link within SU(2) subgroups without changing the action, complementing heatbath rather than serving as the sole update.A trajectory consists of n_hb heatbath sweeps followed by n_ovr overrelaxation sweeps.Source: §2.3 · p. 3

2.4. Checkerboard parallelism and reunitarization

  • Because each staple reads opposite-parity sites, links of one parity in a fixed direction update concurrently; a sweep visits both parities and all D directions.After each trajectory, Gram–Schmidt orthonormalization removes accumulated floating-point drift from link matrices.Source: §2.4 · p. 3

2.5. Observables

  • The code measures the average plaquette and Polyakov loop, accumulating spatial–spatial and spatial–temporal plaquettes separately.The Polyakov loop multiplies links along the temporal direction, serially per spatial site while reducing across sites in parallel.Source: §2.5 · pp. 3–4

3. Software design and Kokkos mapping

kwqft uses compile-time parameters and Kokkos abstractions to express its gauge kernels across execution spaces. Its data layout, distributed-memory layer, and host SIMD path support platform-specific execution without changing the physics kernels.

3.1. Scope and build configuration

  • Compile-time N and D settings parameterize the gauge group and dimension, while configure-time options select Serial, OpenMP, CUDA, HIP, or SYCL execution.Compile-time constants allow loops to be unrolled or vectorized and small matrices to remain in registers, but require a rebuild for each (N, D) pair; SYCL was configured but not benchmarked.Source: §3.1 · p. 4

3.2. Gauge field representation

  • The structure-of-arrays layout makes each matrix element contiguous across links, supporting coalesced GPU access and contiguous CPU SIMD batches.The full format stores all N² complex elements; the SU(3) two-row format stores two rows and reconstructs the third.Source: §3.2 · p. 4

3.3. Indexing, neighbor access, and genericity

  • Templated coordinates and lazy shift views resolve neighbors on access, avoiding neighbor tables and allowing staple address arithmetic to be fused.MPI shift resolution can return a local address or a received ghost-buffer address; checkerboarding requires even extents in specified directions.Source: §3.3 · p. 4

3.4. Node level mapping to Kokkos

  • Kokkos range policies update links, parallel_for handles reunitarization, and parallel_reduce measures observables.Each site owns a random stream indexed by its checkerboard site, with the same sequence across backends; timings and operation or memory counts are diagnostics, not claims against hardware peaks.Source: §3.4 · p. 5, Table 1 · p. 5

3.5. Distributed memory decomposition

  • MPI partitions the lattice into rank-owned subdomains and tracks stale ghost blocks so consumers refresh only data requiring communication.Boundary sites update first; nonblocking halo exchange proceeds while interior sites update, then completes before dependent launches.Source: §3.5 · p. 5

3.6. Host SIMD vectorization

  • Host SIMD batches checkerboard sites across lanes, applying staple, link, matrix, heatbath, and overrelaxation operations across the batch.ISA-specialized de-interleaving loads and stores support ARM SVE/NEON and x86 AVX-512/AVX2; accelerator kernels retain one site per device thread.Source: §3.6 · p. 6

4. Verification

kwqft is checked against an exact two-dimensional solution and published three- and four-dimensional plaquette values. The comparisons show close agreement, with one high-precision three-dimensional discrepancy discussed as inconclusive.

4. Verification

  • Verification ensembles use 20,000 trajectories from a cold start, discard the first 10,000, and analyze the remainder with blocked jackknife errors.The reported comparisons were generated with the Kokkos CUDA backend on one NVIDIA A100; smaller-lattice checks found backend results consistent within statistical uncertainties.Source: §4 · p. 6

4.1. Two dimensions: comparison with the exact solution

  • The two-dimensional reference uses factorized single-plaquette integrals, truncating the q sum at |q|≤30 after convergence to machine precision over the considered β range.The exact solvability follows from independent plaquettes in axial gauge, up to exponentially small finite-volume corrections.Source: §4.1 · p. 6
  • Across β=1–100, kwqft closely matches exact plaquette results for SU(2)–SU(5) and SU(15), including strong- and weak-coupling regions.Figure 1 compares measurements on a 100 × 100 lattice with exact curves and shows residuals in its lower panel.Source: §4.1 · p. 6, Figure 1 · p. 7

4.2. Three and four dimensions: comparison with published data

  • In 21 of 22 three- and four-dimensional cases, kwqft agrees with references within 2.3 combined standard deviations; 16 cases agree within one.The exception is SU(8) at β=100 in three dimensions, differing by 3.2 × 10^-6 (about 4.5 combined standard deviations); the authors do not regard it as evidence of systematic error.Source: §4.2 · p. 7
  • The three-dimensional comparisons use Hietanen and Kurkela values for N=2–8, while four-dimensional comparisons use Bali et al. values for N=2–7 and N=17.Three parameter sets have identical central values; the published results do not specify the algorithm or hardware used to produce them.Source: §4.2 · p. 7

5. Performance methodology

The methodology evaluates a common application source across CPU, GPU, and DCU platforms using trajectory-time measurements. Its experiments cover native-CUDA comparison, host SIMD, distributed scaling, and trajectory time across platforms as N varies.

5.1. Platforms

  • The platform set includes Arm processors with NEON or SVE, an x86 processor with AVX-512, NVIDIA GPUs, and a Hygon DCU using HIP.LineShine LX2 runs bind one MPI rank and its OpenMP threads and memory to each cluster’s locality-matched HBM domain.Source: §5.1 · pp. 7–8

5.2. Measurements and metrics

  • Trajectory time is wall-clock time measured on rank 0 between fences around one trajectory’s update sweeps, averaged over the final 50 or 100 trajectories.Runs contain 50–200 trajectories from a cold start, omitting startup effects such as cache-cold launches and MPI connection establishment.Source: §5.2 · p. 9
  • The experiments compare Kokkos with native CUDA on an A100, SIMD with scalar builds on LX2, strong scaling across 16–512 nodes, and trajectory time across platforms as N varies.The cross-platform workload uses N=2–15 on a 24^3 × 64 lattice at β=6; the strong-scaling workload fixes the global volume at 128^3 × 256.Source: §5.2 · p. 9

6. Performance results

The benchmarks show that kwqft combines near-native CUDA performance on an A100 with SIMD gains and MPI+OpenMP scaling on LineShine. Across several architectures, one source supports SU(N) runs with performance patterns that vary by platform and N.

6.1. Kokkos CUDA backend versus native CUDA

  • 81.5% Kokkos efficiency at L=8 rises to 90–93% for L=32–64 and 99.4% at L=72, as the Kokkos and native CUDA curves nearly coincide.The authors attribute the remaining 7–10% gap at saturated sizes to kernel organization; they caution against over-interpreting the near-unity value at L=72.Source: §6.1 · pp. 9–10, Figure 2 · p. 10

6.2. Host SIMD speedup on LineShine

  • 1.7–2.1× faster for N≤10 and 2.4–3.5× faster for N=11–15, the eight-lane SIMD build outperforms the scalar OpenMP build.The speedup grows with N as dense algebra becomes more prominent and the scalar build faces register and instruction-cache pressure.Source: §6.2 · p. 10, Table 5 · p. 11

6.3. MPI+OpenMP scaling on LineShine LX2

  • 19.5 s to 1.2 s at 512 nodes gives 16.2× speedup for 32× more resources, or 51% strong-scaling efficiency.Super-linear gains through 128 nodes coincide with smaller per-rank fields and greater cache and HBM locality; at 512 nodes, halo exchanges can no longer be hidden behind interior updates.Source: §6.3 · pp. 10–11, Figure 3 · p. 11

6.4. Trajectory time versus 𝑁across architectures

  • For N≥7, the LX2 is fastest; the A100 leads at N≤5, while the DCU is within 1.3–2.6× of the A100 for N≥4 and beats the RTX 5060 throughout.Trajectory time increases with N on every platform, approaching N^3 behavior at large N; the RTX 5060 series ends at N=10 because of its 8 GB memory.Source: Figure 4 · p. 12, §6.4 · p. 12

6.5. Interpretation

  • The same source ran across the tested CPU and accelerator platforms without source-level changes, with execution space, SIMD width, and ISA-specific load/store specialization selected at configure time.The experiments span near-native A100 performance, LineShine host SIMD and scaling, and SU(N) scaling across CPU and accelerator architectures.Source: §6.5 · p. 12

7. Conclusions

kwqft implements Wilson pure-gauge Monte Carlo for SU(N) in arbitrary space-time dimensions using Kokkos, with distributed halo exchange and host SIMD support. Its physics verification and performance results support competitive execution from a single processor to a large distributed system.

7. Conclusions

  • Compile-time parameters express SU(N) and dimensional generality, while Kokkos abstractions decouple algorithms from backend-specific launch and memory interfaces.The halo layer overlaps communication with interior updates, and host SIMD accelerates the same update kernels without altering the CUDA, HIP, or SYCL device path.Source: §7 · pp. 12–13
  • 81%–99% of native CUDA throughput on an NVIDIA A100, 1.7–3.5× OpenMP SIMD acceleration on SVE-512, and 16.2× strong-scaling speedup on LineShine support competitive portability.The OpenMP acceleration is for N=2 to 15; the LineShine result is for a 128^3 × 256 SU(4) lattice scaled from 16 to 512 nodes.Source: §7 · p. 13
  • The implementation reproduces the exact two-dimensional plaquette for SU(2)–SU(15) and published three- and four-dimensional values for N up to 17.Source: §7 · p. 13
  • Extending kwqft to dynamical fermions remains a longer-term objective requiring pseudofermions, portable Krylov solvers, force kernels, reversible symplectic integrators, and Hybrid Monte Carlo.The gauge representation, halo layer, and host SIMD infrastructure are designed for reuse in those components.Source: §7 · p. 13

CRediT authorship contribution statement

Wei Sun is credited with conceptualization, methodology, software, validation, investigation, data curation, visualization, and writing.

CRediT authorship contribution statement

Loading 2610.09766v1…