Scaling Lattice QCD beyond 100 GPUs

R. Babich; M. A. Clark; B. Joó; G. Shi; R. C. Brower; S. Gottlieb

Source-linked AI summary

Scaling Lattice QCD beyond 100 GPUs

R. Babich, M. A. Clark, B. Joó, G. Shi, R. C. Brower, S. Gottlieb

arXiv:1109.2935v1Sep 13, 2011hep-latphysics.comp-ph

TL;DR

GPU methods had advanced LQCD analysis, but gauge-field generation still required capability-scale strong scaling at O(100) GPUs. The paper combines multidimensional parallelization with an additive Schwarz preconditioner and demonstrates scaling beyond 100 GPUs, with tests reaching 256 GPUs for Wilson-clover and improved staggered discretizations.

Problem
Gauge-field generation requires strong scaling into the O(100)-GPU regime, which had not previously been achieved for LQCD.
Method
The paper parallelizes the Dirac operator in multiple dimensions and uses an additive Schwarz domain-decomposed preconditioner with GCR.
Results
The Wilson-clover solver was successfully strong-scaled to over 100 GPUs, while experiments used partitions of up to 256 GPUs and included asqtad solvers.
Takeaways & Limitations
GPU clusters can deliver over 10 teraflops and support lattice-ensemble generation, while multidimensional parallelization enables asqtad solvers at leading-edge lattice volumes.
Takeaways & Limitations
Beyond 32 GPUs, scaling departs from ideal because growing communication overhead makes the implementation communications-bound.

Abstract

from arXiv · show

Over the past five years, graphics processing units (GPUs) have had a transformational effect on numerical lattice quantum chromodynamics (LQCD) calculations in nuclear and particle physics. While GPUs have been applied with great success to the post-Monte Carlo "analysis" phase which accounts for a substantial fraction of the workload in a typical LQCD calculation, the initial Monte Carlo "gauge field generation" phase requires capability-level supercomputing, corresponding to O(100) GPUs or more. Such strong scaling has not been previously achieved. In this contribution, we demonstrate that using a multi-dimensional parallelization strategy and a domain-decomposed preconditioner allows us to scale into this regime. We present results for two popular discretizations of the Dirac operator, Wilson-clover and improved staggered, employing up to 256 GPUs on the Edge cluster at Lawrence Livermore National Laboratory.

1. INTRODUCTION

LQCD GPU codes must overcome communication and strong-scaling limits to support the O(100)-GPU capability computing required for gauge-field generation. This work combines multidimensional parallelization with an additive Schwarz preconditioner and evaluates the approach through 256 GPUs.

GPU LQCD solvers achieve multi-teraflop performance comparable to QCD codes running on capability machines with several thousand cores.
PCI-E and heterogeneous-system imbalances can make closely coupled nearest-neighbor stencil codes communications-bound at large GPU counts.
O(100) GPUs are required for the gauge-generation phase because its Monte Carlo Markov chains demand strong scaling.
Previous implementations parallelized only in time, limiting GPU count and motivating additional dimensions and communication-reducing domain decomposition.
The paper parallelizes QUDA's discretized Dirac operator across multiple dimensions, investigates an additive Schwarz preconditioner for GCR, and tests partitions up to 256 GPUs.

2. LATTICE QCD

Lattice QCD converts non-perturbative QCD calculations into repeated sparse linear solves on discretized spacetime lattices. This section introduces the Wilson-clover and asqtad discretizations and their distinct field and stencil structures.

Lattice QCD provides a model-independent, non-perturbative method for QCD where low-energy perturbation theory fails.
Gauge-field configurations are generated sequentially, while observables are evaluated in a task-parallel analysis phase across ensemble configurations.
Finite-difference discretization turns the Dirac operator into a large sparse matrix, reducing quark calculations to repeated linear systems.
The study focuses on Wilson-clover and asqtad, two widely used discretizations with different treatments of fermion doublers and discretization errors.
Wilson-clover acts on 12-component color-spinors and includes link matrices, spin projectors, a clover correction, and the quark mass.
Asqtad reduces spin degrees of freedom per site to one, while smeared and extended gauge fields reduce its discretization errors.

3. ITERATIVE SOLVERS

LQCD's large sparse systems are solved iteratively, but strong scaling increases communication relative to local computation and eventually makes traditional Krylov methods communication-bound. Additive Schwarz preconditioning addresses this constraint by solving independent subdomains within flexible GCR.

Wilson-clover and staggered formulations require repeated sparse matrix-vector products within iterative Krylov solvers.
Physical quark masses produce nearly indefinite matrices, while linear solvers consume 80–99% of execution time in both calculation stages.
As core counts rise at fixed lattice volume, increasing surface-to-volume ratios raise communication costs until local computation can no longer hide communication.
3.2 Additive Schwarz preconditioner: Additive Schwarz partitions the domain into blocks with Dirichlet boundaries, allowing independent local solves without inter-block communication.
3.2 Additive Schwarz preconditioner: Because local Schwarz systems are solved iteratively, additive Schwarz is used as a preconditioner for the flexible GCR outer solver.

4. OVERVIEW OF RELATED WORK

Earlier LQCD GPU work largely targeted single GPUs or small clusters, whereas this paper targets O(100)-GPU scaling with a QDR Infiniband interconnect. Its approach builds on prior domain-decomposition work while extending the scale substantially.

GPU LQCD research rapidly expanded after CUDA, including multi-GPU work on Wilson, overlap, and domain-wall formulations.
Most reported implementations remained single-GPU, with relatively little literature on multi-GPU LQCD beyond selected QUDA and overlap implementations.
Prior restricted additive Schwarz work used a four-node Gigabit Ethernet cluster, while this paper targets O(100) GPUs over QDR Infiniband.

5. QUDA

QUDA provides optimized CUDA kernels and wrapper code for major LQCD workloads, using memory-access optimizations and supporting multi-GPU applications. Earlier multi-GPU work partitioned only the time dimension, limiting strong scaling for realistic volumes.

QUDA: QUDA supplies optimized CUDA kernels and wrapper code for time-consuming LQCD components, including gauge smearing and force-term computations.It interfaces with the Chroma and MILC applications.
QUDA: Single-GPU QUDA solvers reach up to 24% of peak performance through thread-per-site execution, reordered fields, memory coalescing, and texture-cache use.The strategy assigns each GPU thread responsibility for updating one lattice site.
QUDA: Partitioning only the time dimension enabled strong and weak scaling to 32 GPUs but severely limited strong scaling for realistic lattice volumes.The limitation arose from increasing surface-to-volume ratio as GPU counts grew.

6. MULTI-DIMENSIONAL PARTITIONING

The implementation partitions the lattice across multiple dimensions, exchanges boundary data through ghost zones, and overlaps computation with communication using CUDA streams. As subvolumes shrink, communication increasingly dominates computation, motivating communication-reduction techniques.

6.1 General strategy: Each GPU receives a four-dimensional subvolume bounded by up to eight three-dimensional faces, with boundary data stored in ghost zones.Kernels read either local subvolume data or the corresponding neighbor buffer.
6.1 General strategy: Ghost-zone exchange occurs only along partitioned dimensions, while padding reduces partition camping and also provides gauge-field ghost-zone storage.Gauge-field ghost zones are transferred once at the beginning of a solve.
6.2 Interior and exterior kernels: The Dirac operator uses one exterior kernel for every partitioned dimension, extending the earlier interior–boundary decomposition beyond time-only partitioning.The general case therefore uses up to four exterior kernels.
6.3 Computation, communication, and streams: CUDA streams overlap computation with communication using two streams per dimension plus one stream for interior and exterior kernels, totaling 9 streams.Gather kernels launch immediately so communication can begin in all directions.
6.3 Computation, communication, and streams: For small subvolumes, communication can exceed interior-kernel runtime, leaving the GPU idle and degrading performance.The implementation used extra host memory copies because GPU-pinned memory was incompatible with MPI-pinned memory on the study cluster.

7. DIRAC OPERATOR PERFORMANCE

Performance experiments on the Edge cluster evaluate strong scaling for Wilson-clover and asqtad operators across multiple precisions and partitioning strategies. Scaling becomes communication-bound at high GPU counts, while multidimensional partitioning improves asqtad performance at 256 GPUs.

7.2 Wilson-clover: 256 GPUs were used for the Wilson-clover strong-scaling study, which showed significant departures from ideal scaling beyond 32 GPUs.The study used volume V = 32^3 × 256 with 12 gauge reconstruction.
7.2 Wilson-clover: Beyond 32 GPUs, Wilson-clover becomes communications-bound as increasing surface-to-volume ratio reduces the ability to overlap computation with communication.The performance advantage of half precision over single precision also diminishes as communication overhead grows.
7.3 Improved staggered: For asqtad, faster kernels matter more at low GPU counts, whereas minimizing surface-to-volume ratio becomes increasingly important as GPU count rises.The performance study used volume V = 64^3 × 192 with no gauge reconstruction.
7.3 Improved staggered: At 256 GPUs, the XYZT partitioning scheme achieves the best asqtad performance despite having the worst single-GPU performance.The result reflects the increasing importance of partition geometry at higher GPU counts.

8. BUILDING SCALABLE SOLVERS

The section develops communication-reducing solvers for large GPU counts, combining domain decomposition with mixed precision and restarted GCR. These choices reduce communication and memory pressure while preserving the required solution accuracy.

Domain-decomposed preconditioning: GCR-DD switches off inter-GPU communication within each domain, restricting boundary contributions and reductions to local subdomains.The additive Schwarz preconditioner makes each domain solve independently.
Mixed-precision GCR-DD: The mixed-precision GCR-DD algorithm builds and orthogonalizes a bounded Krylov space before restarting when a size or convergence condition is reached.The maximum Krylov-space size is controlled by kmax, alongside residual-based restart conditions.
Mixed-precision GCR-DD: Mixed precision solves the preconditioned system in half precision, while GCR restarts the Krylov process in high precision.This uses low precision for the Krylov space and high precision for correction and restarting.
Improved-staggered solver strategy: Multi-shift methods cannot be restarted and require storing multiple solution and direction vectors, increasing bandwidth pressure and memory requirements.These restrictions make the mixed-precision strategy less straightforward and motivate sequential refinement.
Improved-staggered solver strategy: Multi-shift CG is implemented in single precision, followed by mixed-precision sequential CG refinement of each solution vector to the desired tolerance.The strategy retains simultaneous multi-shift solving while correcting the resulting solutions sequentially.
Improved-staggered solver strategy: The modified multi-shift strategy achieves double-precision accuracy with most operations in single precision and a shorter time to solution than sequential mixed-precision CG.Its sustained flop rate is lower because of additional linear algebra, but the total solution time is significantly shorter.

9. SOLVER PERFORMANCE RESULTS

The Wilson-clover GCR-DD solver scales to 256 GPUs and reduces time to solution beyond 32 GPUs, while the improved staggered multi-shift solver achieves strong scaling from 64 to 256 GPUs.

Wilson-clover: The Wilson-clover GCR-DD solver scales to 256 GPUs and improves time to solution over BiCGstab by 1.52x, 1.63x, and 1.64x at 64, 128, and 256 GPUs, respectively.At 32 GPUs, BiCGstab remains the superior solver.
Scaling limits: At 256 GPUs, the GCR-DD scaling reaches the algorithm’s limit because a significant fraction of computation still requires full communication.The resulting Amdahl’s law effect produces identical slowdown slopes for GCR and BiCGstab from 128 to 256 GPUs.
Wilson-clover: GCR-DD achieves greater than 10 Tflops on partitions of 128 GPUs and above, comparable to the 10-17 Tflops attained by leadership-class systems on larger core partitions.The comparison uses the same lattice volume and double-precision accuracy for the leadership-system benchmarks.

10. CONCLUSIONS

The paper demonstrates strong scaling of Wilson-clover LQCD beyond 100 GPUs and enables GPU use for asqtad solvers at leading-edge lattice volumes. It identifies the simple non-overlapping Schwarz preconditioner as an initial step toward improved scaling on heterogeneous architectures.

Conclusions: Multi-dimensional parallelization and an additive Schwarz preconditioner strong-scale the Wilson-clover solver to over 100 GPUs.The resulting GPU clusters deliver in excess of 10 teraflops, described as the minimum capability needed for lattice-ensemble generation.
Conclusions: Multi-dimensional parallelization enables GPUs to run asqtad solvers at leading-edge lattice volumes, previously prevented by the operator’s decreased locality.
Future improvements: The simple non-overlapping additive Schwarz preconditioner is only a first step, with overlapping domains or multilevel blocking proposed to improve scaling substantially.These alternatives are intended to exploit multiple levels of memory locality in GPU clusters.
Future improvements: The authors frame GPUs and Schwarz preconditioning as components of a broader restructuring for heterogeneous architectures with deep memory hierarchies.They anticipate domain decomposition, mixed-precision solvers, and data compression among the relevant tools for exascale LQCD.