PaperlayerView exact revision on arXiv

Source-linked AI summary

Solving Lattice QCD systems of equations using mixed precision solvers on GPUs

M. A. Clark, R. Babich, K. Barros, R. C. Brower, C. Rebbi

arXiv:0911.3191v2hep-lat

TL;DR

Lattice QCD requires efficient solution of large sparse Dirac systems, whose matrix-vector products dominate computation. The paper implements Wilson-Dirac kernels and reliable-update mixed-precision Krylov solvers on NVIDIA CUDA, finding that half-precision BiCGstab with double-precision reliable updates gives the least time to solution and solver performance above 100 Gflops.

  • Problem

    Lattice QCD requires efficient solution of many large sparse linear systems involving the Dirac operator, with matrix-vector application as the most compute-intensive kernel.

  • Method

    The paper implements optimized Wilson-Dirac matrix-vector routines and mixed-precision Krylov solvers using reliable updates on NVIDIA CUDA.

  • Results

    The half-precision BiCGstab solver with double-precision reliable updates gives the least time to solution, while solver performance exceeds 100 Gflops.

  • Takeaways & Limitations

    Reliable-update mixed precision provides full double-precision accuracy while using single or half precision for the bulk of computation, outperforming defect-correction in the tested systems.

  • Takeaways & Limitations

    The reported work is an initial single-GPU exploration, with multi-GPU adaptation left for future work.

Abstract

from arXiv · show

Modern graphics hardware is designed for highly parallel numerical tasks and promises significant cost and performance benefits for many scientific applications. One such application is lattice quantum chromodyamics (lattice QCD), where the main computational challenge is to efficiently solve the discretized Dirac equation in the presence of an SU(3) gauge field. Using NVIDIA's CUDA platform we have implemented a Wilson-Dirac sparse matrix-vector product that performs at up to 40 Gflops, 135 Gflops and 212 Gflops for double, single and half precision respectively on NVIDIA's GeForce GTX 280 GPU. We have developed a new mixed precision approach for Krylov solvers using reliable updates which allows for full double precision accuracy while using only single or half precision arithmetic for the bulk of the computation. The resulting BiCGstab and CG solvers run in excess of 100 Gflops and, in terms of iterations until convergence, perform better than the usual defect-correction approach for mixed precision.

1. Introduction

GPUs offer highly parallel computation and high memory bandwidth for scientific workloads, motivating their use for lattice QCD. The paper targets GPU implementation of lattice-QCD solvers and introduces reliable-update mixed precision to retain double precision accuracy without explicit restarts.

  • Motivation: GPUs emphasize parallelism and memory bandwidth, offering a different trade-off from commodity CPUs for highly parallel scientific workloads.The paper frames GPUs as a potentially better balance for intrinsically parallel tasks than processors optimized for clock rate, caches, and instruction-level parallelism.
  • Motivation: Lattice QCD requires many linear solves involving the large sparse Dirac operator, making efficient solvers a grand-challenge problem.Applying the Dirac operator to a vector is the most compute-intensive kernel in these solvers.
  • Approach: The implementation uses NVIDIA CUDA and the C for CUDA language to implement lattice-QCD computations on GPUs.CUDA provides direct access to the NVIDIA GPU through a C-like programming model.
  • Mixed precision: GPU performance depends strongly on precision because single precision can be much faster than double precision, especially on GPUs.The reported hardware difference can reach an order of magnitude when double precision is supported.
  • Mixed precision: Reliable updates replace the iterated residual periodically with a high-precision true residual and apply high-precision groupwise solution updates without explicit restarts.The method is designed to preserve the Krylov search space while using lower precision for most computation.

2. Graphics Processing Units

The paper describes NVIDIA GPU hardware and CUDA’s thread-based programming model, emphasizing occupancy, memory behavior, and synchronization constraints. These characteristics motivate keeping computation on the GPU and carefully managing memory and precision.

  • Hardware: The GTX 280 contains 30 multiprocessors and 240 cores, with latency hidden through high multiprocessor occupancy.Device memory is shared across multiprocessors and has relatively high latency.
  • Hardware: GPU shared memory is manually managed, much faster than device memory, so minimizing device-memory access is important.Unlike typical CPUs, GPUs do not contain a large memory cache.
  • Data movement: The entire solver is implemented on the GPU because PCI Express transfers can impose a severe bandwidth constraint.The initial matrix and source download and final solution upload are negligible compared with solver operation time.
  • CUDA programming model: CUDA launches large numbers of parallel threads organized into blocks and grids, with shared memory available within each block.A lattice-QCD application can assign one thread to each lattice site.
  • CUDA programming model: CUDA permits synchronization within a thread block but not between different blocks during one kernel invocation.This restriction shapes how kernels and inter-thread communication must be organized.
  • CUDA programming model: High occupancy reduces latency penalties, while divergent execution paths serialize the paths across a 32-thread warp.The GPU dynamically schedules thread blocks, and groups of 32 threads are intended to follow the same execution path.

3. The Wilson-Dirac Matrix

The Wilson-Dirac matrix models lattice-QCD quark propagation as a large sparse linear system over spacetime, spin, and color degrees of freedom. The paper uses Krylov solvers and exploits even-odd structure while focusing GPU optimization on the Dslash matrix-vector operation.

  • Matrix formulation: The Wilson-Dirac matrix discretizes the Dirac operator and adds a scaled Laplace term to remove spurious fermion doublers.It acts on a tensor product of four-dimensional Euclidean spacetime, spin space, and color space.
  • Matrix formulation: Each spacetime point contributes a complex-valued 12-component color-spinor, and the complete vector is called a spinor field.The gauge field consists of SU(3) link matrices acting in color space.
  • Linear systems: Quark physics requires many solutions for different source vectors and gauge fields, so the paper considers Krylov sub-space iterative solvers.The matrix dimensions and sparsity make iterative methods the feasible choice.
  • Even-odd preconditioning: Even-odd ordering produces a Schur complement whose condition number is 2-3 times smaller than that of the full matrix.Solving the Schur complement also halves vector-memory requirements, although applying it requires an extra temporary field when preserving the original source.
  • GPU kernel: The Dslash operation is the application of D to a spinor field and constitutes the most compute-intensive part of a Krylov solver.Efficient GPU solving therefore depends on implementing this matrix-vector product efficiently.
  • GPU kernel: Generic sparse-matrix GPU libraries cannot exploit Wilson-Dirac symmetries, limiting their performance relative to structure-aware implementations.The cited library reaches up to 30 Gflops versus a 933 Gflops GTX 280 peak, while Wilson-Dirac symmetries can reduce memory throughput.

4. CUDA Implementation

The CUDA implementation specializes the Wilson-Dirac Dslash kernel around lattice structure, memory access, gauge-field representations, and precision-specific constraints. These choices reduce bandwidth demands and achieve up to 135 Gflops in single precision, 40 Gflops in double precision, and 212 Gflops in pseudo half precision on a GTX 280.

  • Kernel and data layout: The Dslash operation requires 1,440 bytes per site, making its performance severely limited by memory bandwidth.Each site requires 1,320 floating-point operations, 8(24 + 18) float loads, and 24 float saves.
  • Kernel and data layout: Gauge fixing, γ4 diagonalization, and 8-number gauge parametrization reduce Dslash traffic from 1,440 to 960 bytes.The combined strategy reduces bandwidth requirements by one third.
  • Single precision: 135 Gflops is the peak performance of the 8-parameter single-precision kernel with temporal gauge fixing.The 8-parameter reconstruction adds operations, but its overhead is not noticeable because the kernel is bandwidth starved.
  • Double precision: 40 Gflops is the peak of the gauge-fixed 12-parameter double-precision kernel, nearly three times below the equivalent single-precision kernel.Double precision also doubles memory traffic and increases register pressure, limiting occupancy to 128 active threads.
  • Half precision: 212 Gflops is the peak of the 12-parameter pseudo half-precision kernel, 50% faster than the fastest single-precision kernel.For pseudo half precision, packing, unpacking, and gauge reconstruction double the matrix-vector product cost, leaving performance floating-point limited.
  • Performance versus accuracy: Across precisions, the 12-parameter reconstruction is at least one to two orders of magnitude more numerically accurate than the 8-parameter reconstruction.The 8-parameter method performs more operations, including trigonometric functions and division by a subtracted quantity.

5. Mixed Precision Krylov Solvers

The paper implements CG and BiCGstab solvers across precisions on GPUs and compares defect-correction with reliable updates for mixed-precision convergence. Reliable updates preserve precision without restarting, generally improving iteration behavior and time to solution, though memory use and light-mass behavior shape the choice.

  • 5.1. Implementation and Performance: CG and BiCGstab are implemented for double, single, and half precision, with separate CUDA kernels for linear-algebra operations and matrix-vector products.CG uses the normal equations because the Wilson-Dirac matrix is not Hermitian positive definite; BiCGstab applies directly.
  • 5.2. Defect-Correction: Defect-correction performs low-precision inner solves with high-precision residual calculation and solution accumulation, but discards the Krylov subspace after each restart.This discarded search space can increase the total iterations required for convergence.
  • 5.3. Reliable Updates: Reliable updates replace the iterated residual with the true residual when it decreases by a factor δ, maintaining precision without the overhead of restarting.The method uses low-precision updates and periodically applies high-precision residual corrections; δ must satisfy ulp_in < δ < 1.
  • 5.4. Results: At light quark masses, reliable-update single precision increases iterations by at most 15%, while half precision converges at the lightest mass with a 34% increase for δ = 10^-1.Larger δ values generally produce iteration counts closer to the double-precision baseline.
  • 5.4. Results: Defect-correction can be preferable for memory-constrained problems, whereas reliable updates are preferred at light quark masses when memory is not a concern.Reliable updates require additional low-precision spinor fields and gauge-field storage for both precisions.
  • 5.4. Results: The reliable-update Single 8 solver delivers an almost constant factor 3 speedup, while Half 12 achieves factors of 4 to 3.5 as quark mass decreases.Relative to the approximately 39 Gflops double-precision solver, these results correspond to effective solver performance well in excess of 100 Gflops.

6. Conclusions

The paper implements optimized Wilson-Dirac matrix-vector routines and GPU Krylov solvers, emphasizing reduced memory bandwidth. Reliable-update mixed precision solvers achieve the best tested time to solution, while the resulting GPU performance offers cost and power advantages.

  • Conclusions: Optimized Wilson-Dirac matrix-vector routines form the basis of CUDA Krylov solvers, with high performance obtained by minimizing memory bandwidth.The approach accepts increased floating-point operation counts to reduce data movement.
  • Conclusions: The half precision BiCGstab solver with double precision reliable updates gave the least time to solution for the tested systems.The method was compared against conventional defect-correction mixed precision.
  • Conclusions: At 100 Gflops sustained inverter performance, GPUs offer considerable cost and power savings over traditional cluster and massively parallel architectures.The authors expect GPUs to become more attractive for QCD calculations as new solutions appear.
  • Conclusions: The solver became the mainstay of the open-source QUDA library and was interfaced with Chroma, CPS, and QDP/C.The library also supports the Wilson-Clover discretization in addition to the Wilson-Dirac matrix.
  • Conclusions: Figure 7 compares time to solution and speedup versus double precision for BiCGstab with δ = 0.1, ǫ = 10^-12, and volume = 243 × 64.Solid lines represent time to solution, while dashed lines represent speedup.
  • Conclusions: Future work includes implementing additional costly lattice-QCD operations and adapting the code for multi-GPU environments.The multi-GPU performance was left for future reporting.

Appendix A. SU(3) Matrix Reconstruction

The appendix introduces notation for representing a general SU(3) matrix by naming its matrix components. The components are arranged as three rows of complex-valued entries.

  • Appendix A. SU(3) Matrix Reconstruction: A general SU(3) matrix is introduced by labeling its components with the symbols a1 through a3, b1 through b3, and c1 through c3.These labels provide the notation used for subsequent reconstruction steps.
  • Appendix A. SU(3) Matrix Reconstruction: The notation identifies nine matrix entries before the appendix describes reduced representations and reconstruction procedures.The supplied passage establishes the component labels but does not specify the reduction itself.
  • Appendix A. SU(3) Matrix Reconstruction: The labeled components are grouped into three rows, beginning with a1 a2 a3 and followed by b1 b2 b3 and c1 c2 c3.The row-wise labeling supports references to stored and reconstructed rows.

Appendix A.1. 12 number parametrization

The 12-number parametrization reconstructs an SU(3) matrix after storing only its first two rows. The third row is then obtained from those stored rows.

  • Appendix A.1. 12 number parametrization: Only the first two rows of the SU(3) matrix are stored in the 12-number parametrization.The passage presents this as the starting point for reconstructing the remaining row.
  • Appendix A.1. 12 number parametrization: The third row is reconstructed from the stored first two rows rather than explicitly stored.This reduces the data that must be retained for the full matrix.
  • Appendix A.1. 12 number parametrization: The parametrization therefore represents the full SU(3) matrix through its first two stored rows and a derived third row.The supplied passage states the storage and reconstruction relationship without giving the reconstruction formula.

Appendix A.2. 8 number parametrization

The 8-number parametrization reduces SU(3) link storage by reconstructing matrix information from selected components and phases. It is designed for efficient, numerically stable GPU evaluation, with a remaining singularity at |a1| = 1.

  • Appendix A.2. 8 number parametrization: The modified reconstruction reduces the matrix representation to eight real numbers after intermediate parametrization and normalization constraints.The method stores phases of a1 and c1 rather than using stereographic projection.
  • Appendix A.2. 8 number parametrization: The row vectors b and c lie in the plane orthogonal to a, allowing the original SU(3) matrix to be obtained from an SU(2) rotation of that plane.The reconstruction uses the relation between the row vectors and the orthogonal plane.
  • Appendix A.2. 8 number parametrization: Storing the phases of a1 and c1 avoids the two singularities introduced by the earlier stereographic-projection approach.The full numbers are recovered through trigonometric functions.
  • Appendix A.2. 8 number parametrization: The phase-based reconstruction is numerically stable and efficient on GPUs because it uses fast trigonometric and square-root functions.The passage explicitly connects the stability and efficiency to the reconstruction design and GPU hardware.
  • Appendix A.2. 8 number parametrization: In half precision, setting a negative square-root argument to zero substantially improves numerical stability.This is a half-precision-specific safeguard.
  • Appendix A.2. 8 number parametrization: The reconstruction still has a singularity at |a1| = 1 because of the normalization factor N.This is the stated remaining limitation of the method.

Appendix B. Dirac Matrix Conventions

The appendix contrasts chiral and non-relativistic spin-projector bases for Wilson-Dirac matrix operations. The UKQCD basis reduces temporal-gather memory traffic, while transformations preserve compatibility with lattice QCD packages.

  • The DeGrand-Rossi basis is a conventional chiral basis for the spin projectors in Wilson-Dirac off-diagonals.
  • In the chiral basis, applying a projector requires loading all spinor components regardless of dimension or direction.
  • The UKQCD basis requires only upper or lower spin components for temporal backward or forward gathers, respectively.
  • This selective loading halves temporal-gather memory traffic and increases kernel performance.
  • Basis transformations are applied during host-to-device transfer and reversed during device-to-host transfer for compatibility with existing lattice QCD packages.
Loading 0911.3191v2…