Calculation of HELAS amplitudes for QCD processes using graphics processing unit (GPU)

We use a graphics processing unit (GPU) for fast calculations of helicity amplitudes of quark and gluon scattering processes in massless QCD. New HEGET ({\bf H}ELAS {\bf E}valuation with {\bf G}PU {\bf E}nhanced {\bf T}echnology) codes for gluon self-interactions are introduced, and a C++ program to convert the MadGraph generated FORTRAN codes into HEGET codes in CUDA (a C-platform for general purpose computing on GPU) is created. Because of the proliferation of the number of Feynman diagrams and the number of independent color amplitudes, the maximum number of final state jets we can evaluate on a GPU is limited to 4 for pure gluon processes ($gg\to 4g$), or 5 for processes with one or more quark lines such as $q\bar{q}\to 5g$ and $qq\to qq+3g$. Compared with the usual CPU-based programs, we obtain 60-100 times better performance on the GPU, except for 5-jet production processes and the $gg\to 4g$ processes for which the GPU gain over the CPU is about 20.


Introduction
In our previous report [1] we introduced a C-language [2] version of the HELAS codes [3], HEGET (HELAS Evaluation with GPU Enhanced Technology), which can be used to compute helicity amplitudes on a GPU (Graphics Processing Unit). Encouraging results with 40-150 times faster computation speed over the CPU performance were obtained for pure QED processes, qq → nγ, for n = 2 to 8 in pp collisions.
In this paper, we extend our study to QCD processes with massless quarks and gluons. The HEGET routines for massless quarks and gluons are identical to those of quarks and photons introduced in [1], and the qqg vertex function structure is also the same as the qqγ functions. The only new additional routines are those for ggg and gggg vertices. For the QED processes studied in ref. [1], we found that the present CUDA compiler cannot process qq → 6γ amplitude with 6! ≈ 700 Feynman diagrams, and we need to subdivide the HEGET codes into small pieces for 6γ and 7γ processes. In the case of 8γ production with 8! ≈ 4 × 10 4 Feynman diagrams, we have not been able to compile the program even after subdivision into small pieces. We also encountered serious slow down when the program accesses global memory during the parallel processing period. Therefore, our concern for evaluating the a e-mail: junichi.kanzaki@kek.jp b e-mail: naotoshi@post.kek.jp c e-mail: tstelzer@uiuc.edu QCD processes on a GPU is the proliferation of the number of diagrams, as well as the number of independent color amplitudes which come with different color weights.
The paper is organized as follows. In section 2, we present the cross section formula for n-jet production processes in pp collisions in the quark-parton model, or in the leading order of perturbative QCD with scale-dependent parton distribution functions (PDF's). In section 3, we review briefly the structure of GPU computing by using HEGET codes, and give basic parameters of the GPU and CPU machines used in this analysis. In section 4, we introduce new HEGET functions for ggg and gggg vertices. Section 5 gives our results and section 6 summarizes our findings. Appendix lists all the new HEGET codes introduced in section 4.
2 Physics Process 2.1 n-jet production in pp collisions The cross section for n-jet production processes can be expressed as dσ = {a,b} dx a dx b D a/p (x a ,Q) D b/p (x b ,Q) dσ(ŝ) , (1) where D a/p and D b/p are the scale (Q) dependent parton distribution functions (PDF's), x a and x b are the momentum fractions of the partons a and b, respectively, in the Table 1. The number of Feynman diagrams and the color bases for QCD processes studied in this paper.
No. of jets gg → gluons uu → gluons uu → uu+gluons in the final state #diagrams #colors #diagrams #colors #diagrams #colors  2  6  6  3  2  2  2  3  45  24  18  6  10  8  4  510  120  159  24  76  40  5  7245  720  1890  120  786  240 right-and left-moving protons. For the total pp collision energy of gives the invariant mass squared of the hard collision process The subprocess cross section is computed in the leading order as where is the invariant n-body phase space, λ i are the helicities of the initial and final partons, n a and n b are the color degree of freedom of the initial partons, a and b, respectively, and c i represents the color indices of the initial and final partons. When there are more than one gluons or identical quarks in the final states, an appropriate statistical factor should be multiplied on the phase space dΦ n in eq. (5). The Helicity amplitudes for the process (1) can be expressed as where the summation is over all the Feynman diagrams. The subscripts λ i stand for a given combination of helicities (±1 for both quarks and gluons in the HELAS convention [3]), and the subscripts c i correspond to a set of color indices (1, 2, 3 for flowing-IN quarks, 1, 2, 3 for flowing-OUT quarks, and 1 to 8 for gluons). In MadGraph [4] the amplitudes are expanded as in the color bases T ci α which are made from the SU(3) generators in the fundamental representation [5] The color factors are computed as where n a,b = 3 for q and q, n a,b = 8 for gluons, and the summation is over all {c i} = {c a , c b , c 1 , . . . , c n }. The color sum-averaged square amplitudes are computed as The cross sections are then expressed as where we introduce the helicity sum-average symbol as In this paper the following three types of multi-jet production processes are computed: The number of contributing Feynman diagrams and the number of color bases for the above processes are summarized in Table 1, which includes those for the process, gg → 5g. We note here that the number of diagrams (7245) for gg → 5g exceeds that of the uu → 7γ process (7! ≈ 5040), for which we could run the converted MadGraph codes on a GPU, only after division into small pieces [1]. In fact, we have not been able to run the gg → 5g program on GPU even after dividing the program into more than 100 pieces; as explained in section 5.4. Proliferation of the number of independent color basis vectors is also a serious concern for GPU computing, since the color matrix N of eq. (9) has m(m + 1)/2 elements when there are m independent basis vectors T ci a . For example, the process uu → uuggg has m = 240 color basis vectors from Table 1, and the matrix has 3×10 4 elements. The matrix exceeding 16000 elements cannot be stored in the 64kB constant memory, while storing it in the global memory will result in serious loss of efficiency in parallel computing. Therefore, the method to handle summation over color degrees of freedom is a serious concern in GPU computing.

Selection criteria for jets
Total and differential cross sections of the processes (13) in pp collisions at √ s = 14TeV are computed in this paper. We introduce final state cuts for all the jets as follows: where η i and p Ti are the rapidity and the transverse momentum of the i-th jet, respectively, in the pp collisions rest frame along the right-moving (p z = |p|) proton momentum direction, and p Tij is the relative transverse momentum [6] between the jets i and j defined by Here ∆R ij measures the boost-invariant angular separation between the jets. As for the parton distribution function (PDF), we use the set CTEQ6L1 [7] and the factorization scale is chosen to be the cut-off p T value, Q = p cut T = 20 GeV. The QCD coupling constant is also fixed as which is obtained from the MS coupling at Q = m Z , α s (m Z ) MS = 0.118 [8] by using the NLO renormalization group equations with 5-flavors.
3 Computation on the GPU

GPU and its host PC
For the computation of the cross sections of QCD n-jet production processes we use the same GPU and host PC as in the previous report [1]. In particular we use a GeForce GTX280 by NVIDIA [11] with 240 processors, whose parameters are summarized in Table 2. It is controlled by a Linux PC with Fedora 8 on a CPU whose properties are summarized in Table 3. Programs which are used for the computation of the cross sections are developed with the CUDA [2] environment introduced by NVIDIA [11] for general purpose GPU computing.

Program structure
Our program computes the total cross sections and distributions of the QCD n-jet production processes via the following procedure: 1. initialization of the program, 2. random number generation for multiple phase-space points {p a , p b , p 1 , . . . , p n } and helicities {λ i } on the CPU, 3. transfer of random numbers to the GPU, 4. generation of helicities and momenta of initial and final partons using random numbers, and compute amplitudes (J λi ) a of eq. (8) for all the color bases on the GPU, 5. multiplying the amplitudes and their complex conjugate with the color matrix N ab of eq. (9) and summing them up as in eq. (10), and multiply the PDF's of the incoming partons on the GPU, 6. transferring momenta and helicities for external particles, computed weights and the color summed squared amplitudes to the CPU, and 7. summing up all values to obtain the total cross section and distributions on the CPU.
Program steps between the generation of random numbers (2) and the summation of computed cross sections (7) are repeated until we obtain sufficient statistics for the cross section and all distributions.

Color matrix calculation
In order to compute the cross sections of the QCD multijet production processes, multiplications of the large color matrix N ab of eq. (9), the vector of color-bases amplitudes (J λi ) α of eq. (8) and its complex conjugate have to be performed, as in eq. (10). For large n-jet processes, like gg → 4 gluons, uu → 5 gluons and uu → uu+3 gluons, the dimensions of color matrices exceed 100, and the number of multiplication becomes larger than 10 4 . These matrices cannot be stored in the constant memory (64kB for the GTX280; see Table 2) which is accessed in parallel, while storing them in the global memory (1GB for GTX280) results in serious slow-down of the GPU. We find that multiplications for the color-summation in eq. (10) can be reduced significantly as follows. The color matrix of eq. (9) contains many elements with the same value. We count the number of different non-zero elements in the color matrix and find the results shown in Table 4. We find for instance that among the 240×(240+1)/2 = 28, 720 elements of the color matrix for the uu → uu + 3g process, there are only 60 unique ones.
In general, the number of different elements in the color matrix grows linearly rather than quadratically as the number of color basis vectors grows. Since the numbers in Table 4 are small enough, we can store them in the constant memory which is accessed quickly by each parallel processor.
Before we arrive at the above solution adopted in this study, we examined the possibility of summing over colors via Monte Carlo. Let us briefly report, in passing, on this exercise.
In the Monte Carlo color summation approach, we evaluate the matrix element M ci λi (7) for a given set of momenta {p i }, helicities {λ i } and colors {c i }, and sum the squared amplitudes over randomly generated sets of {p i , λ i , c i }. This method turns out not to be efficient because in the color basis using the fundamental representation of the SU(3) generators adopted by MadGraph, most of the basis vectors T ci a vanish for a given color configuration {c i }. As an example, gg → 4g has 5! = 120 color basis vectors (see Table 1), which take the form for the configuration {c i } = (a 1 , a 2 , . . . , a 6 ) where a i denotes the color index of the gluon i taking an integer value between 1 and 8. Among the 8 6 ≈ 260, 000 configurations, only 12% give non-zero values. Moreover, as many as 75% of the color configurations give vanishing results for all the 120 basis vectors. Although the efficiency can be improved by changing the color basis, we find that our solution of evaluating the exact summation over colors is superior to the Monte Carlo summation method for all the processes which we report in this paper.

New HEGET functions
The HEGET functions for massless quarks and gluons are the same as those introduced in the previous report [1].
The qqg vertex functions are identical to the qqγ functions of ref.
[1] except for the coupling constant; for the vertex where g s = √ 4πα s is the strong coupling constant and T a ij is an SU (3)) generator in the fundamental representation. For example, the qqg vertex function is computed by the HEGET function iovxx0 as iovxx0(cmplx* fi, cmplx* fo, cmplx* vc, float g, cmplx vertex vertex vertex) (20) where the coupling constants are following the convention of MadGraph [4] and the color amplitude In the rest of this section, we introduce new HEGET functions for three-vector boson (VVV) and four-vector boson (VVVV) vertices. All the new HEGET functions are listed in Table 5, and their contents are given in Appendix. Also shown in Table 5 is the correspondence between the HEGET functions and the HELAS subroutines [3].

VVV: three vector boson vertex
For the ggg vertex we introduce two HEGET functions, vvvxxx and jvvxx0. They correspond to HELAS subroutines, VVVXXX, and JVVXXX, respectively, for massless particles; see Table 5.

vvvxxx
The HEGET function vvvxxx (List 1 in Appendix) computes the amplitude of the VVV vertex from vector boson wave functions, whether they are on-shell or off-shell. The function has the arguments: vvvxxx(cmplx* ga, cmplx* gb, cmplx* gc, float g, cmplx vertex vertex vertex) where the inputs and the outputs are: Inputs: cmplx ga [6] wavefunction of gluon with color index, a cmplx gb [6] wavefunction of gluon with color index, b cmplx gc [6]   The coupling constant is in the HEGET function (24), following the convention of MadGraph [4]. In order to reproduce the amplitudes associated with the ggg vertex Lagrangian of eq. (23), the color factor associated with the ggg vertex is if abc . More explicitly, the vertex amplitude for eq. (23) is 1 by using the output (vertex vertex vertex) in eq. (24). Also note the HELAS convention [3] of using the flowing-OUT momenta and quantum numbers for all bosons.

jvvxx0
This HEGET function jvvxx0 (List 2 in Appendix) computes the off-shell vector wavefunction from the threepoint gauge boson coupling in eq. (23). The vector propagator is given in the Feynman gauge for a massless vector bosons like gluons. It has the arguments: jvvxx0(cmplx* ga, cmplx* gb, float g, cmplx* jvv jvv jvv) where the inputs and the outputs are: Inputs: cmplx ga [6] wavefunction of gluon with color index, a cmplx gb [6] wavefunction of gluon with color index, b float g coupling constant of the VVV vertex Outputs: cmplx jvv [6] vector current j µ (gc : ga, gb) which has a color index, c (29) As in eq. (27) the color amplitude for the off-shell current is if abc (jvv jvv jvv) .

VVVV: four vector boson vertex
For the ggggg vertex we introduce two HEGET functions, ggggxx and jgggx0, listed in Table 5. They correspond to HELAS subroutines, GGGGXX and JGGGXX, respectively, for massless particles.

ggggxx
The HEGET function ggggxx (List 3 in Appendix) computes the portion of the amplitude of the gggg amplitude where the first and the third, and hence also, the second and the fourth gluon wave functions are contracted, whether the gluons are on-shell or off-shell. The function has the arguments: ggggxx(cmplx* ga, cmplx* gb, cmplx* gc, cmplx* gd, float g, cmplx vertex) where the inputs and the outputs are: Inputs: cmplx ga [6] wavefunction of gluon with color index, a cmplx gb [6] wavefunction of gluon with color index, b cmplx gc [6] wavefunction of gluon with color index, c cmplx gd [6] wavefunction of gluon with color index, d float gg coupling constant of VVV vertex Outputs: cmplx vertex amplitude of the VVVV vertex The coupling constant gg for the gggg vertex is In order to obtain the complete amplitude, the function must be called three times (once for each color structure) with the following permutations: ggggxx(ga, gb, gc, gd, gg, v1 v1 v1) (35a) ggggxx(ga, gc, gd, gb, gg, v2 v2 v2) (35b) ggggxx(ga, gd, gb, gc, gg, v3 v3 v3) The color amplitudes are then expressed as

jgggx0
The HEGET function jgggx0 (List 4 in Appendix) computes an off-shell gluon current from the four-point gluon coupling, including the gluon propagator in the Feynman gauge. It has the arguments: jgggx0(cmplx* ga, cmplx* gb, cmplx* gc, float gg, cmplx* jggg) where the inputs and the outputs are: Inputs: cmplx ga [6] wavefunction of gluon with color index, a cmplx gb [6] wavefunction of gluon with color index, b cmplx gc [6] wavefunction of gluon with color index, c float gg coupling constants of the VVVV vertex Outputs: cmplx jggg [6] vector current j µ (gd : ga, gb, gc) which has a color index, d (38) The function (37) computes off-shell gluon wave function with three specific color index d which comes along with a specific color factor. As in eq. (35) it should be called three times jgggx0(ga,gb,gc,gg,j1 j1 j1) (39a) jgggx0(gc,ga,gb,gg,j2 j2 j2) (39b) jgggx0(gb,gc,ga,gg,j3 j3 j3) (39c) to give the off-shell gluon with the color factor

Comparison of total cross sections
In order to validate the new HEGET functions which are introduced in this report, we compare the total cross sections of n-jet production processes computed on the GPU with those calculated by other programs which are based on the FORTRAN version of the HELAS library. We use MadGraph/MadEvent [4] and another independent FOR-TRAN program which uses the Monte Carlo integration program, BASES [12], as references. Due to the limited support for the double precision computation capabilities on the GPU, the whole computations with HEGET on a GTX280 are done with single precision, while the other programs with HELAS in FORTRAN compute cross sections with double precision. For the calculation of the n-jet production cross sections we use the same physics parameters as the Mad-Graph/MadEvent for all programs, and the same final state cuts of eq. (14) for all processes and all programs. The parton distribution functions of CTEQ6L1 [7] and the same factorization and renormalization scales, Q = p cut T = 20GeV, are also used.
Results for the computation of the total cross sections are summarized in Tables 6, 7 and 8 for gg → gluons, uu → gluons and uu → uu+gluons, respectively. We find the results obtained by the HEGET functions agree with those from the other programs within the statistics of generated number of events.
We note that multi-jet events that satisfy the final state cuts of eq. (14), where all jets are in the central region in |η| < 2.5 (14a) and their transverse momentum about the beam direction (14b) and among each other (14c) greater than 20 GeV, are dominated by pure gluonic processes in Table 6. The cross sections for uu → ng process in Table 7 are small because of uu annihilation. We note that the crossing-related non-annihilation processes, ug → u + (n − 1)g, have exactly the same number of diagrams and color bases, hence can be evaluated with essentially the same amount of computation time.

Comparison of the processing time
As already described in our previous report [1], we prepare two versions of the programs in the same structure for the computation of the total cross sections. One is written in CUDA, a C-based language, and can be executed on the GPU. The other is written in C and can be executed on the CPU. Using a standard C library function we measure the time between the start of the transfer of random numbers to the GPU and the end of the transfer of computed results back to the CPU.
In Fig. 1, the measured process time in µsec for one event of n-jet production processes is shown for the GPU (GTX280) and the CPU (Linux PC with Fedora 8). They are plotted against the number of jets in the final state. Because the process time per event on the GPU depends [1] strongly on the number of allocated registers at the compilation by the CUDA and the size of thread blocks at the execution time, we scan combination of these parameters for the fastest event process time on the GPU.
The upper three lines in Fig. 1 show the event process times on the CPU. They correspond to gg → n-jets denoted as gg, uu → n-jets as uu and uu → uu + n-jets as uu, respectively. For processes with small numbers of jets, e.g. n jet = 2, the event process times for different processes are all around 4.5 µsec. This is probably because they are dominated by computation steps other than the amplitude calculations, such as computations of the PDF factors and the data transfer between GPU and CPU, which are common to all physics processes. When the number of jets becomes larger, the event process time for the same number jets in the final states is roughly proportional to the number of diagrams of each process listed in Table 1.
The lower three lines in Fig. 1 show the event process times on a GTX280. They also correspond to gg → n-jets denoted as gg, uu → n-jets as uu and uu → uu + n-jets as uu, respectively. As the number of jets becomes larger, the process time on the GPU grows more rapidly than that on the CPU. For the n jet = 4 case, the event process time of gg → 4 gluons is larger than the expected time from the proportionality to the number of diagrams of the other processes, uu → 4 gluons and uu → uu+2 gluons. In other words, the event process time on GPU grows faster than what we expect from the growth of the number of Feynman diagrams.
For instance, the event process times ratio for gg → 4g and gg → 3g on the CPU are roughly 120 µsec/14 µsec  ∼ 8.6, which roughly agrees with the ratio of the numbers of Feynman diagrams (Table 1), 510/45 ∼ 11. The corresponding ratio on GPU is 3.8 µsec/0.1 µsec ∼ 38, which is significantly larger.
For the same number of jets, we also observe that the event process times on the CPU are roughly proportional to the number of diagrams. Fir n jet = 4, the ratio of the process times for gg → 4g to uu → 4g are about 120 µsec/29 µsec ∼ 4.1 on CPU, as compared to the ratio of the number of Feynman diagrams in Table 1, 510/159 ∼ 3.2. The same applies to n jet = 5 between uu → 5g and uu → uuggg, where Feynman diagrams have the ratio 1890/786 ∼ 2.4 from Table 1, and the event process time on the CPU gives 300 µsec/180 µsec ∼ 1.7, also in rough agreement.
On the other hand, the event process times on the GPU for gg → 4g and uu → 4g have a ratio 3.8 µsec/0.45 µsec ∼ 8.4 which is much larger than the ratio of the diagram numbers; while that for uu → 5g and uu → uuggg has the ratio of 11µsec/9.5µsec ∼ 1.15. Although we do not fully understand the above behavior of the event process time on the GPU, we find that they tends to scale as the product of the number of Feynman diagrams and the number of color bases, while the event process times on the CPU are not sensitive to the latter. This is probably because as the number of color bases grows, more amplitudes, (J λi ) α in eq. (8), should be stored and then called to compute the color sum, eq. (10). These observations tell us that the relative weight of the color matrix computation in the GPU computing is very significant even after identifying the independent elements of the color matrix N αβ in eq. (9) as listed in Table 4.

Comparison of performance of GPU and CPU
The ratios of event process times between CPU and GPU are shown in Fig. 2. Three lines correspond to gg → n-jets denoted as gg, uu → n-jets as uu and uu → uu+(n−2)jets as uu, respectively. The performance ratios exceed 100 for the processes with small numbers of jets (n jet ≤ 3) in the final state. For n jet = 4 and 5, the performance ratios gradually drop to less than 40. For processes with large numbers of color bases, the ratios are smaller. For gg → 4 gluons, which has 120 color bases, the ratio is about 30, and for uu → uu+3 gluons, which has 240 color bases, the ratio becomes about 20.

Note on gg → 5g study
Among five-jet production processes we have not been able to run the program for gg → 5g. This process has 7245 di- agrams and 720 color basis vectors. In order to compile the program for the computation of this process, we use the technique developed in the previous study [1]. By dividing the program into about 140 pieces we were able to compile the gg → 5g program. Compilation takes about 90 min. on a Linux PC. The total size of the compiled program exceeds 200 MB, and we were not able to execute this compiled program on a GTX280.

Summary
We have shown the results of our attempt to evaluate QCD multi-jet production processes at hadron colliders on a GPU [11], Graphic Processing Unit, following the encouraging results obtained for QED multi-photon production processes in ref. [1].
Our achievements and findings may be summarized as follows.
-A new set of HEGET functions written in CUDA [2], a C-language platform developed by NVIDIA for general purpose GPU computing, are introduced to compute triple and quartic gluon vertices. The HEGET routines for massless quarks were introduced in ref. [1], and the routine for photons [1] can be used for gluons. In addition, the HEGET functions for the qqg vertex are the same as those for the qqγ vertex introduced in ref. [1]. -The HELAS amplitude code generated by MadGraph [4] is converted to a CUDA program which calls HEGET functions for the following three type of subprocesses: gg → ng (n ≤ 5), uu → ng (n ≤ 5), and uu → uu + ng (n ≤ 3). -Summation over color degrees of freedom was performed on a GPU by identifying the same valued elements of the color matrix of eq. (9), in order to reduce the memory size. -All the HEGET programs for up to 5 jets passed the CUDA compiler after division into small pieces. However, we could not execute the program for the process gg → 5g. Accordingly, comparisons of performance be-tween GPU and CPU are done for the multi-jet production processes up to 5 jets, excluding the purely gluonic subprocess. -Event process times of the GPU program on GTX280 are more than 100 times faster than the CPU program for all the processes up to 3-jets, while the gain is reduced to 60 for 4-jets with one or two quark lines, and to 30 for the purely gluonic process. It further goes down to 30 and 20 for 5-jet production processes with one and two quark lines, respectively. -We find that one cause of the rapid loss of GPU gain over CPU as the number of jets increases is the growth in the number of color bases. GPU programs slow down for processes with larger numbers of color basis vectors, while the performance of the CPU programs is not affected much. -All computations on the GPU were performed with single precision accuracy. A factor of 2.5 to 4 slower performance is found for double precision computation on the GPU.  [4].re); p1[1] = (ga [5].re); p1[2] = (ga [5].im); p1[3] = (ga [4].im); p2[0] = (gb [4].re); p2[1] = (gb [5].re); p2[2] = (gb [5].im); p2[3] = (gb [4].im); q[0] = -(jvv [4].re); q[1] = -(jvv [5].re); q[2] = -(jvv [5].im); q[3] = -(jvv [4].im);