Baryonic and mesonic 3-point functions with open spin indices

We have implemented a new way of computing three-point correlation functions. It is based on a factorization of the entire correlation function into two parts which are evaluated with open spin- (and to some extent flavor-) indices. This allows us to estimate the two contributions simultaneously for many different initial and final states and momenta, with little computational overhead. We explain this factorization as well as its efficient implementation in a new library which has been written to provide the necessary functionality on modern parallel architectures and on CPUs, including Intel's Xeon Phi series.


Introduction
Observables constructed with the use of three-point correlation functions can describe a multitude of physical phenomena, such as the parton structure of hadrons or their weak transitions, depending on the operator type at the insertion. The relevant strong interaction matrix elements can be computed using Lattice Quantum Chromodynamics. Traditionally, a method employed to that aim was the sequential source method [1]. The numerical cost of this is high because a new inversion is necessary for each final momentum at the sink. In this contribution we present a stochastic algorithm which circumvents this limitation and disentangles the number of inversions of the lattice Dirac operator from the number of sink/insertion momenta. The implementation we propose parallelizes the computations in such a way that multiple source positions and multiple insertion positions can be estimated simultaneously. Moreover, by storing the uncontracted data, with all spin indices open, on disk, we enable the user to analyse any channel of interest at a later stage.
An aspect of our implementation which we wish to highlight in this contribution is the extensive use of vectorization. As the compute power of today's processors relies on longer and longer vector registers of additional vector processing units, an adequate data layout is required to efficiently use these resources. We explain our data layout which supports large vector registers, therefore enabling us to efficiently run our code on the Intel's Xeon Phi processors featuring 512 bit AVX vector instructions.
Our implementation is based on a framework developed specifically for Intel's Xeon Phi processors, called LibHadronAnalysis see, e.g., [2,3]. This library contains additional routines for computing meson and baryon spectra, meson and baryon distribution amplitudes and many other objects. Compared to a naive implementation in Chroma [4] using QDP++ [4] objects it provides speed-up factors of the order 10-20 on the KNC and KNL architectures.
Below we describe in detail the algorithm. In particular we explain how the computation of a threepoint correlation function can be separated into two, largely independent parts, the "spectator" and the "insertion" parts. We pay special attention to the parallelization schemes which are different for each of these parts. Subsequently, we present benchmarks showing the performance of our implementation on a KNL cluster, and then we conclude.

Stochastic baryonic and mesonic three-point correlation functions
In this section we introduce the three-point correlation function we are interested in and show how it can be factorized into two, largely independent parts. We only show explicit formulae for the case of meson three-point functions, for the sake of notational simplicity. A similar approach was used in Refs. [5][6][7][8], however, here we keep all spin indices open.

General structure
A three-point meson correlation function (c.f., figure 1) with a source C(r), a sink A(x ) and an insertion operator I(y), located at timeslices r 4 , x 4 and y 4 respectively, reads are the annihilation, creation and insertion operators respectively, with f i ∈ {l, s, c} (light, strange, charm). G f i (r, x) is a standard fermion propagator from x to r of flavor f i . We use the convention to denote the annihilation spin and color operator indices with primed Greek and Latin letters α , a , creation operator indices with ordinary Greek and Latin letters α, a, and the insertion operator indices with tilde Greek and Latin lettersα,ã. Γ ins can contain local derivatives and A and C may contain quark smearing. At this point we replace one of the propagators by its stochastic estimate where the sum runs over N realizations of the noise source vector η i (x ), with the properties The (η i ) (x) are time partitioned and set to zero, unless In figure 1 we show the generic structure of a three-point correlation function in the meson case. The central ellipse denotes a meson created with operator Eq. (3) in the middle of the lattice at timeslice r 4 . The meson is annihilated by the operator from Eq. (2) at the timeslice x 4 as depicted by the left-most ellipse. Arrows denote exact point-to-all propagators. The star at timeslice y 4 denotes one of the possible positions of the insertion operator. The wiggly line is used to plot the stochastic all-to-all propagator. We call these four elements the "forward" correlation function as opposed to the shaded mirror-reflected graph on the right hand side which corresponds to the "backward" process. The forward and backward diagrams are estimated simultaneously which allows for increased statistics.
The introduction of the stochastic propagator allows us to factorize the correlation function A(x ) I(y) C(r) into two parts [5] as follows We define the spectator where we have assumed r = 0. Otherwise we have to replace x → x − r, y → y − r.

Spectator part
The computation of the spectator part consists of the contractions of propagators at the timeslices where the source and the sink are located. Naively only the MPI ranks working on the timeslices r 4 and x 4 , x 4 would work. In our implementation we prepare a set of propagators sourced from different temporal source positions, as shown on figure 2. We redistribute these propagators among the different MPI ranks in such a way that each rank has at least one propagator. The computation of the spectator part for each source position is than performed simultaneously. The Fourier transformation Eq. (9) fixes the momentum p at the sink.

Insertion part
The insertion part corresponds to the contraction of the stochastic propagator, i.e., the solution of the lattice Dirac equation sourced by random noise vectors, with the point-to-all propagator. This has to be repeated for each position y 4 of the insertion operator between the sinks x 4 and x 4 and the source r 4 . Again, in order to keep a good workload balance we redistribute the data from the timeslices where the insertion is present, denoted with stars on figure 3, among all MPI ranks in such a way that each rank has approximately the same number of insertion positions to work with, and we perform all the computations in parallel. Note that a separate Fourier transformation in Eq. (10) allows to select a desired momentum q flowing through the insertion.
x 4 τ r 4 x 4 · · · y 44 · · · y 49 · · · Figure 3. Parallelization of the insertion part of the three-point correlation function. Data from timeslices denoted with yellow stars is redistributed among all MPI ranks, such that all ranks have a similar workload.

Performance
We paid special attention to the data layout to enable the use of vectorization. The crutial point was to reorder the many loops in the algorithm. We show this explicitly for pseudocode in listings 1 and 2. The LibHadronAnalysis library [2] was incorporated into the QDP++/Chroma software stack [4], together with the multigrid solver [9][10][11][12][13] optimized for the KNL architecture.
The computations were carried out on a single configuration of the CLS ensemble H101 [14]. Running on 4 − 32 KNLs and distributing the 256 hardware threads on each KNL to 8 MPI tasks yields the strong scaling for the meson spectator and insertion part contraction as shown in figure 4. Note that the computation is done for 50 independently seeded stochastic estimators in forward and backward direction, i.e., the mean values are averaged over 100 computations each. The creation times for the noise and the solution vectors are not considered.
For the spectator/insertion a minimal computation time is achieved using 8/32 KNLs with 8 tasks per node. In this setup the Intel Omni-Path connection between the nodes becomes saturated, the computation is well parallelized and the overhead due to internal communication is not dominant.
In addition the wallclock-time for a particular measurement using two different ranges of momenta on one H101 configuration is evaluated -again 8 KNLs with 8 tasks per node are used. The timings for k = k 2 = 0 and for k , k 2 = 0, . . . , 8 for spectator and insertion momenta are shown in figure 5, where p i = k i 2π/L, q i = k i 2π/L where the integers k i and k i label the momentum components of p and q within the Fourier transformation. Note that in these timings also a baryon measurement is included. In both cases the overall computation time is almost the same LibHadronAnalysis wallclock-time for k 2 = k 2 = 0: ≈ 530s LibHadronAnalysis wallclock-time for k 2 , k 2 = 0, . . . , 8 (93 · 93 mom. combinations): ≈ 575s.
Hence it is possible to produce data for a large number of final momenta without increasing the computation time significantly. In addition the data-layout presented in Eq. (9) and Eq. (10) provides analysis capabilities for various physical channels since the Γ-structures of the source and sink interpolators are not specified during the simulations.
The k 2 , k 2 = 0, . . . , 8 measurement was also performed using the sequential source method [1] which yields a run-time of ≈ 930s. Compared to the above wallclock-time of ≈ 575s this is a speed-up of ≈ 1.6 where we have not taken into account that the stochastic code gives 16 · 16 Γ-combinations at source and sink for free. First tests have shown that the computation time of the analysis code needed to obtain final three-point function results is in the range of a few seconds and therefore is negligible.
Collectively this means that one needs ≈ 82 KNL core hours to perform the above measurement on a single configuration of the H101 assuming that 8 nodes with 64 cores per node are used. Altogether 2016 configurations are available to analyse the H101, i.e., at most ≈ 165000 KNL core hours are needed to analyze the entire ensemble for every combination of source and sink meson and baryon interpolators on a single source time slice and for the given source sink separation ∆t = 10 a. Due to the distribution shown in figure 5 the overall computation time should remain almost constant when increasing the number of source positions, always provided that the computation of the spectator part can be fully parallelized.

Conclusions
The timings presented in the last section reveal that our implementation overtakes currently implemented methods at least by a factor of 1.5 in terms of computation time for high momentum square for a given source-sink structure. However, exploiting the great flexibility of LibHadronAnalysis output data makes it feasible to reuse the data for many source-sink Dirac Γ-matrix combinations which results in an incredible economization. Furthermore it is sufficient to compute the insertion part only once and use it in both, baryon and meson, measurements. Since it is now possible to generate data containing an enormous amount of information it is also necessary to process the data further to finally get physical observables. The analysis software package is still in development and will be released as soon as possible.