Fast and flexible implementations of Wilson, Brillouin and Susskind fermions in lattice QCD

TL;DR

This paper presents a modern, flexible Fortran implementation of Wilson, Brillouin, and Susskind fermion operators in lattice QCD, optimized for modern many-core architectures with SIMD and parallelization.

Contribution

It introduces a portable, high-performance Fortran code applying multiple Dirac operators simultaneously, with various memory layouts, for efficient lattice QCD computations.

Findings

Achieves high performance on modern architectures (up to 880 Gflop/s)

Demonstrates the importance of memory layout optimization for KNL architecture

Provides a self-contained code with Krylov solvers for lattice QCD applications

Abstract

A modern Fortran implementation of three Dirac operators (Wilson, Brillouin, Susskind) in lattice QCD is presented, based on OpenMP shared-memory parallelization and SIMD pragmas. The main idea is to apply a Dirac operator to $N_v$ vectors simultaneously, to ease the memory bandwidth bottleneck. All index computations are left to the compiler and maximum weight is given to portability and flexibility. The lattice volume, $N_x N_y N_z N_t$, the number of colors, $N_c$, and the number of right-hand sides, $N_v$, are parameters defined at compile time. Several memory layout options are compared. The code performs well on modern many-core architectures (480 Gflop/s, 880 Gflop/s, and 780 Gflop/s with $N_v=12$ for the three operators in single precision on a 72-core KNL processor, a $2\times24$-core Skylake node yields similar results). Explicit run-time tests with CG/BiCGstab inverters confirm that the memory layout is relevant for the KNL, but less so for the Skylake architecture. The ancillary code distribution contains all routines, including the single, double, and mixed precision Krylov space solvers, to render it self-contained and ready-to-use.