Lessons learnt from vectorizing Madgraph5_aMC@NLO and porting it to GPUs

A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 1 / 26 Lessons learnt from vectorizing Madgraph5_aMC@NLO and porting it to GPUs (NB! This talk focuses on computing issues - performance i.e. generation speed! - not physics) Andrea Valassi (CERN IT) Thanks to the whole madgraph4gpu development team! Thanks to S. Frixione, Z. Was and the Les Houches 2023 participants for many useful discussions! ECFA Higgs Factories: 2nd Topical Meeting on Generators –Bruxelles, 22nd June 2023 https://indico.cern.ch/event/1266492 [An update over the previous talk at the ECFA 1st Topical Meeting in November 2021 - https://indico.cern.ch/event/1078675] A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 2 / 26 Outline •Motivation and overview •Some results –Performance: throughout speedups on CPU SIMD and on GPUs for LO processes –Functionality: development status, usability for the experiments •Some lessons learnt –Applicability to other (existing and future) Monte Carlo generators –Do’s and dont’s •Some future prospects and challenges –QCD (EWK?) NLO and beyond –Non-technical challenges •Conclusions A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 3 / 26 Motivation 1: Monte Carlo Event Generators in WLCG computing •HL-LHC computing needs are predicted to outpace resource growth: need R&D to improve software •MC generators are essential for HL-LHC physics and use ~5-20% of ATLAS/CMS WLCG CPU budgets –Speeding them up helps address the resource gap and may allow more complex (N)NLO multi-jet simulations •Challenges and opportunities to improve MC software have been discussed in the HSF generator WG –See the WG review paper prepared for the LHCC review in 2020: https://doi.org/10.1007/s41781-021-00055-1 WLCG meeting with LHCC referees, Feb. 2020 This plot is probably obsolete by now! Sherpa speedups >>2 were reported at ACAT A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 4 / 26 Motivation 2: GPUs and vector CPUs are underexploited in HEP •GPUs provide most of the compute power in recent HPCs (e.g. Summit: 95%) –Supercomputers at HPC centers are already heavily used by the LHC experiments on an “opportunistic” basis –But only a small share of HEP software workloads can run on GPUs today –NB: at Higgs factories, maybe speed will not be an issue? But the need to use GPUs will probably be there! •Most WLCG CPUs support wide vector registers (SSE4.2, AVX2 or above) –But only a small share of HEP software workloads efficiently exploit CPU vectorization today •These architectures are best suited to lockstep processing with limited branching... as in MC generators! –MC generators are ideal candidates to exploit data parallelism in GPUs (SIMT) and vector CPUs (SIMD)! –In this talk I describe how we achieved this for Madgraph5_aMC@NLO (via event-level data parallelism) Summit: NVidia V100 GPUs Aurora: Intel Xe GPUs Leonardo: NVidia A100 GPUs Juwels: NVidia A100 GPUs LUMI: AMD MI250X GPUs A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 5 / 26 Madgraph5_aMC@NLO (MG5aMC) •One of the workhorses for event generation in ATLAS and CMS! •MG5aMC production version is in Fortran –Software outer shell: Madevent (random sampling, integration and event generation + I/O, multi-jet merging...) –Software inner core: Matrix Element (ME) calculation code, automatically generated for each physics process •Matrix Element calculations take 95%+ of the CPU time for complex processes (e.g. gg→tҧ tggg ) •And ME calculations are precisely one component that can be “easily” accelerated on GPUs and vector CPUs... https://doi.org/10.1007/JHEP07(2014)079 A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 6 / 26 MG5aMC and the madgraph4gpu project https://doi.org/10.1051/epjconf/202125103045 https://doi.org/10.22323/1.414.0212 PoS(ICHEP2022)212 https://arxiv.org/abs/2303.18244 •madgraph4gpu: speed up Matrix Element calculation in MG5aMC on GPUs and vector CPUs –Collaboration of theoretical/experimental physicists and software engineers –born in HSF generator WG –Extensive details may be found in the vCHEP2021, ICHEP2022 and ACAT2022 conference proceedings •See also the (Feb 2023) Computing Accelerator Forum presentation, https://indico.cern.ch/event/1207838 •Two parallel approaches to reimplement the ME calculation –(1) “CUDACPP”, initial CUDA/C++ targeting NVidia GPUs and SIMD on CPUs (now extended to AMD GPUs using HIP) –(2) SYCL (and other Portability Framework: Kokkos, Alpaka) later addition supporting many GPUs/CPUs (now extended to SIMD) •Several successive development steps over time (in each of cudacpp and PFs) –(a) standalone applications, our initial prototype –we still use this to optimize the ME calculation alone –(b) MadEvent-integrated applications, ~final goal –existing Fortran framework, with faster ME! –(c) orchestration of many MadEvent applications –usable by LHC experiments, same interface with faster ME! A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 7 / 26 MG5aMC: old and new architecture designs 1. STANDALONE (TOY APPLICATIONS) MULTI-EVENT API 2. NEW MADEVENT (GOAL: LHC PROD) MULTI-EVENT API OLD MADEVENT (NOW: LHC PROD) SINGLE-EVENT API First we developed the new ME engines in standalone applications (Amdahl...) SCALAR: NEW BOTTLENECK? PARALLEL: MUCH FASTER! MATRIX ELEMENT: CPU BOTTLENECK IN OLD MADEVENT MATRIX ELEMENTS CUDA/C++ or PFs: cuRAND CUDA/C++ or PFs: RAMBO CUDA/C++ or PFs: MEKERNELS MOMENTA FORTRAN: RANMAR FORTRAN: MADEVENT CUDA/C++ or PFs: MEKERNELS MOMENTA MATRIX ELEMENTS FORTRAN: RANMAR FORTRAN: MADEVENT FORTRAN: MATRIX1 MOMENTA MATRIX ELEMENTS Then we modified the existing all-Fortran MadEvent into a multi-event framework and we injected the new MEs into it Compute Accelerator Forum, February 2023 https://indico.cern.ch/event/1207838 A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 8 / 26 MG5AMC+cudacpp: CUDA/C++, Fortran, bash, python... 1. STANDALONE TOY APPLICATION OK! (2020-2021) 2. MADEVENT (ONE APPLICATION) OK! (2022) 3. MADEVENT (N x APPLICATIONS) ./bin/generate_events BEING TESTED (June 2023) MG5AMCNLO GITHUB + MADGRAPH4GPU GITHUB MG5AMCNLO GITHUB 4. COMPLETE WORKFLOW INCLUDING CODE GENERATION generate.. output.. launch SOON! (2023) A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 9 / 26 ACAT2022 MadEvent with vectorized C++ for gg→tҧ tgg (on a single CPU core) ME speedup ~ x8 (double) and x16 (float) over scalar Fortran Our ME engine reaches the maximum theoretical SIMD speedup! Overall speedup so far~ x6 (double) and x10 (float) over scalar Fortran (Amdahl’s law) 512y = AVX512, ymm registers 512z = AVX512, zmm registers The latter is only better on nodes with 2 FMA units (here an Intel Gold 6148) Compute Accelerator Forum, February 2023 https://indico.cern.ch/event/1207838 A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 16 / 26 •(1) Design computational units with well-defined inputs and outputs! –Beware of hidden inputs and outputs from common blocks and static data... •(2) Keep data parallelism in mind from the start: move from single-event APIs to multi-event APIs! –Well-defined input array of many events, well-defined output array of many events If you design a new Monte Carlo from scratch, these are MUST's, not SHOULD's! Do's and dont's - two simple recommendations PROCESS ONE EVENT 1 IN 1 OUT PROCESS N EVENTS N IN N OUT REENTRANT FUNCTION (NO STATE! THREAD SAFE!) IN OUT STATEFUL FUNCTION "IN" ? "OUT" ? /COMMON/... /COMMON/... static ... static ... An additional technicality: prefer Structure-of-Array (SOA) memory layouts for the inputs and outputs! [Strictly needed only internally for SIMD and useful for GPUs, but good to have also in the API of the function] A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 17 / 26 •You will still need to loop over multiple sets of N events –And the internal implementation of N-event processing may still involve some loops! •N should be at least as big as the number needed to "fill" the GPU or the CPU vector register –On a vector CPU: at most 16 floats in a 512-bit AVX512 register –On a GPU: we need ~16k threads for complex processes (500k for simpler ones!) •NB: I focus on event-level parallelism, but other options exist –In MG5AMC we will investigate using 1 GPU thread per helicity per event... What about loops? And how many are N events? PROCESS N EVENTS N IN N OUT (N x M EVENTS) LOOP OVER 1...M "Process N events": three implementation examples (there can be more!) 1. CPU scalar: internally loop over N events, process each one individually 2. CPU vector: hold the events in a SIMD vector of size N, 3. GPU kernel: each of the N events is processed by one of N GPU threads A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 18 / 26 Helicity amplitudes –same code in CUDA and in vectorized C++ •Old slide! The new code is different, the idea is the same! •Formally the same code for CUDA and scalar/vector C++ –hide type behind a typedef –add a few missing operators SIMD in CUDA/C++ uses compiler vector extensions! Flexible design: being reused also for vectorized SYCL! Automatically generated! A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 19 / 26 MG5AMC is not alone –SHERPA on GPU (BlockGen) •Note: unlike MG5aMC, based on Feynman diagrams, SHERPA uses ~Berends-Giele recursion relations –Allows computations with more final-state jets •No ongoing effort on CPU vectorization (yet) •Planned Les Houches project: a detailed comparison of software performances of MG5AMC and SHERPA –Tentative process list: pp to tt(0-3jets) or Z(0-3jets) –Previously, an old wish of the HSF generator WG –(NB: not a comparison of physics results or distributions) From http://dx.doi.org/10.21468/SciPostPhysCodeb.3 More recent results were presented this week in Les Houches by Max Knobbe A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 20 / 26 NLO, loops Marco Zaro –https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Pavia2015 B, V, R: matrix elements MC: parton shower Sand Hevents: two separate sets of events (different matrix elements) Integral = S+H is positive –but individual events can have negative weights MC@NLO: https://doi.org/10.1088/1126-6708/2002/06/029 Matching NLO QCD and parton showers (avoid double counting) •So far we have only worked on LO QCD processes! •NLO QCD processes are more computationally intensive –More Feynman diagrams –And especially, loop diagrams! (quad precision needed?) –A matching procedure (MC@NLO) must also be applied •We should be able to compute Born and Real emission contributions in our vectorized C++ and CUDA –We should also be able to handle NLO matching using the current MadEvent based infrastructure –The main challenge will be understanding the computational impact of loops (Amdahl bottleneck?) •News (for me!) from some discussions last week at Les Houches –Branching should not be an issue at NLO, but will be at NNLO? Local subtraction schemes... •What the code does depends on where you are in phase space... –NLO and NNLO needs “complicated” functions like polylogarithms (are these supported in SIMD and CUDA?) –Libraries exist to emulate quad precision (even for SIMD and CUDA), we can look at these (strip them down?) •What about EWK beyond-LO corrections? –If I understand correctly, our approach would be portable, and the same types of challenges would apply? Z. Wettersten (+ OM, SR, AV, R. Schoefbeck) A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 21 / 26 Reweighting •Advantages of reweighting: savings in computing costs (no detector simulation), fewer statistical fluctuations •In practice for MG5AMC: read in an LHE file, add weights, write back the modified LHE file –Will use the new matrix element engine in CUDA/C++ – For further details and a status report: Zenny’s poster at CHEP 2023 •Theoretical and technical challenges –NLO reweighting (see O. Mattelaer, https://arxiv.org/abs/1607.00763) –Coverage of phase space in the new parameter set –Reweighting for a given event-by-event helicity and color Z. Wettersten (+ OM, SR, AV, R. Schoefbeck) 1. Generate signal sample at  ref, with wi(  ref)=1 (By definition, background does not depend on ) 2. Full detector simulation (MC truth event properties xi(true) →observed event properties xi) 3. Reweight each event by matrix element ratio Old technique, renewed interest! A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 22 / 26 Reweighting and weight derivatives in parameter estimation •Weight derivative: event-by-event sensitivity to the measured parameter •First: makes it possible to determine the limit error with an ideal detector, and how much (0 to 1) we do worse –with a given luminosity at a FCC-ee, what is the best theoretically achievable measurement on Higgs couplings? • Second: can be used as a basis for an “improved optimal observable” ML method AV Knowing one’s limits: maximum achievable information with an ideal detector -Ideal acceptance, select all signal events Ssel=Stot - Ideal resolution, measured iis that from MC truth (implies ideal rejection of background events,  i=0) Weight Derivative Regression i(MC truth) ~ q( xi(MC) ) Data observable event properties xi(DATA) Fit WDR regressor qi(DATA) = q( xi(DATA) ) qi(MC) = q( xi(MC) ) https://doi.org/10.1051/epjconf/202024506038 https://zenodo.org/record/3715951 A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 23 / 26 Photo by Geran de Klerk on Unsplash A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 24 / 26 Scientific Computing and Software Collaborations (or: working on the bridge between different units and communities) •A big lesson learnt from porting MG5AMC to GPUs: you need collaborations with a mix of skills! •Developing Monte Carlo generator software: which kind of job is this? In which box should it be? –A scientist’s job? A theorist’s job? An experimentalist’s job? A computing engineer’s job? –Do we need dedicated Scientific Computing units in our labs and universities? –Do we need to have dedicated career paths similar to Research Software Engineers? •The challenge: attracting, training, retaining people with the right competencies and interests –Can we attract and motivate young theorists to work on software and computing optimizations? • A theorist colleague I was recently talking to: “We had an opening for working on software optimizations for our Monte Carlo generator. The only suitable candidates were two theorists. But they were concerned that working on software optimizations would harm their future careers as theorists and refused the job. In the end, we did not hire anyone.” –Can we attract and motivate young software engineers to work with us instead of tech or finance companies? I am only reporting a problem here... I do not have a magic-wand solution  A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 25 / 26 ... and finally... •Upcoming: Workshop on software acceleration of MC event generators –Where? at CERN –When? in ~October-November 2023 (any time constraints with important MC events?) –Organised together with LPCC, MCnet, HSF... •Contact me if you are interested and/or if you have any suggestions A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 32 / 26 Q1 Q2 Q3 Q4 2024 Q1 Q2 Q3 Q4 2023 2022 2021 2020 20192018 (3) Outlook: WIP, plans, ideas for more speed and features Q3 Q4 Sep 2022 HEPscore benchmarking workshop Sep 2022 Lugano hackathon (sub-event-level parallelism...) Q1-Q2 alpha release?... From Q3 2022: Zenny’s PhD work (reweighting, NLO...) A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 33 / 26 CUDA/C++: a single source code approach (so far...) •The main difference between our CPU (C++) and GPU (CUDA) implementations is the following –on the CPU, all computations and all memory access takes place on the host –on the GPU, it is necessary to distinguish computations and memory accesses on the host and on the device •Within the GPU code, the amount of code that is specific to NVidia/CUDA is minimal –Memory allocations (cudaMalloc), encapsulated within host/device buffer classes –Kernel executions (<<<...>>>), encapsulated within very few specific classes –A few specific types or features (thrust::complex, curand, cuBLAS), also encapsulated in specific classes •The rest of the code is (at least formally –see example later) ~identical for C++ and CUDA! –There are almost more differences between scalar and vector C++ code... •Therefore, we presently use a single source code approach for CUDA/C++ (with #ifdef __CUDA__) –We might review this later on –as it sometimes imposes slightly unnatural choices, and may hinder readability –But so far it has allowed us to make rapid progress for both CUDA and C++ in parallel! A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 34 / 26 Memory layouts –AOS, SOA, AOSOA We have experimented with three possible memory layouts for momenta (1) Array-of-Structures AOS: momenta[Nevt][Npar][4] (2) Structure-of-Arrays SOA: momenta[Npar][4][Nevt] (3) AOSOA: momenta[Npag][Npar][4][Nepp] with Nevt = Npag (“pages”) * Nepp (“events per page”) We are using AOSOA’s as the current default – but this is still largely configurable •For CPU vectorization, AOSOAs (or SOAs) are absolutely mandatory! –We use an AOSOA with Nepp equal to the SIMD vector size NeppV –and an aligned malloc is needed too! –For performance comparison we also build a no-SIMD mode with Nepp=1, which is effectively an AOS •For GPUs (1 event per thread), AOSOAs are faster (fewer memory accesses) but not strictly necessary –We use Nepp=4(8) for doubles(floats) so that each page is 32 bytes (the “sector” size, or L2 cache line size) –For a given number of “requests”, AOS uses 4 times more “sectors” (transactions) than AOSOA with Nepp=4 •Coding for SIMD is more complex than coding for GPUs... MATRIX ELEMENTS CUDA/C++: MEKERNELS MOMENTA Matrix element calculation (simplified example) –inputs[4*Npar*Nevt] =(x,y,z,E)-momentum of Npar particles for Nevt events (n-dim array, substructure) –outputs[Nevt] = matrix element for Nevt events (1-dim array, no substructure) Example: Npar=6 particles for the 2→4process gg→tҧ tgg A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 35 / 26 Monitoring GPU memory access –NSight Compute •Explicitly collect two relevant profiler metrics in NSight Compute –“requests” : l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum –“sectors” (i.e. transactions, network roundtrips): l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum –this is from old tests in August 2020 (issue #16), the profiler metrics names may have changed since then •Profile AOS against the AOSOA baseline –same number of “requests” in AOS and AOSOA –AOS needs 4 times as many “sectors” as AOSOA (which fits 4 doubles in a 32-byte cache line) –in other words: AOSOA provides coalesced memory access, AOS does not –for what it is worth (not much!), the actual slowdown in this e+e−→+−example was only 7% however A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 36 / 26 In practice in MG5aMC: use helicity amplitudes and QCD color decomposition 1. (for each helicity ) compute partial amplitudes Jffor each color ordering permutation f (sum diagrams relevant to f) 2. (for each helicity ) compute the sum over colors as the quadratic form JCJ* using the constant color matrix C 3. sum over helicities [Example for 𝑔𝑔 → 𝑡ҧ 𝑡𝑔𝑔𝑔: 128 helicities (before and after filtering)] Each step computes many events 𝒑in parallel! CPU: 1 SIMD event-vector at a time. GPU: 1 event per thread. Inside the ME calculation: Feynman diagrams, colors, helicities Given the momenta Ԧ𝑝 of initial+final partons in one specific event Sum over all helicity combinations of initial+final partons Sum over all color combinations c of initial+final partons Include all Feynman diagrams d allowed for the given and c Example for 𝑔𝑔 → 𝑡ҧ 𝑡𝑔𝑔𝑔: 1240 Feynman diagrams (using helicity amplitudes) This takes ~40% of the CPU time for this process Example for 𝑔𝑔 → 𝑡ҧ 𝑡𝑔𝑔𝑔: 120 color ordering permutations, 120x120 matrix This takes ~60% of the CPU time for this process A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 37 / 26 C++ vectorization –why choose Compiler Vector Extensions? •Portable –available in gcc, clang, icpx (from clang) with minimal differences –Do not require any external libraries or tools (VC, VCL, VecCore, xSIMD, UME::SIMD, or SYCL...) •Powerful, but easy to use –No need to debug auto-vectorization when it does not vectorize –As powerful as intrinsics, but much easier to write (higher-level abstractions) •Intuitive –CVEs force you to think in terms of vector types! •Minor disadvantage –no vector complex type out of the box –But it was easy to write it in our case (RRRRIIII memory layout) as we only need + -   –A few extensions for Boolean vector masks were needed, too •One technical detail: we malloc a standard (aligned!) fptype* and reinterpret_cast as fptype_v*... HUGE THANKS TO SEBASTIEN PONCE for his Practical Vectorization lectures mentioning CVEs! A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 38 / 26 Monitoring lockstep –GPU NSight compute, CPU disassemble •GPU: explicitly collect one profiler metric in NSight Compute –“branch efficiency” : sm__sass_average_branch_targets_threads_uniform.pct –old test (May 2021 issue #25) comparing two code bases: no-divergence baseline has 100% efficiency, alternative with minor forced divergence has 96% efficiency (and is 20% slower) •CPU: the best lockstep metric IMO is the speedup over a no-SIMD case (reach theoretical maximum!) –but is also useful to disassemble the object using objdump and categorize SIMD intrinsics symbols... # Symbols in .o SSE4.2 (xmm) AVX2 (ymm) AVX512 (ymm) AVX512 (zmm) Build type Scalar 4534 0 0 0 SSE4.2 12916 0 0 0 AVX2 0 10630 0 0 256-bit AVX512 0 10366 12 0 512-bit AVX512 0 1267 60 9910 4a90ec2 gg→tҧ tgg ACAT2022 A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 39 / 26 Code generation: how did we bootstrap the project? A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 40 / 26 Code generation: from many “epochs” to a single evolving “epoch” Code generation infrastructure -Python framework and “cudacpp” plugin - Fortran, C++, CUDA templates -Post-generation patches (temporary...) Automatically generated code - Fortran framework (Madevent) - CUDA/C++ Matrix Elements (1) develop on top of auto-generated code (2) backport immediately to code generation infrastructure (3) re-generate NEW MODEL (since end 2021) OLD MODEL (2020early 2021) (1) MG5AMC Python framework, Fortran templates: “upstream” https://github.com/mg5amcnlo/mg5amcnlo (2) CUDACPP plugin, post-generation patches, generated CUDA/C++ physics processes: our https://github.com/madgraph5/madgraph4gpu WIP to-do before a release: full port from madgraph4gpu to mg5amcnlo (remove postgeneration Fortran patches, add CUDACPP upstream) ... and beyond A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 41 / 26 Why focus on complex processes? Compute >> memory! •We are lucky: the more complex the physics process, the less relevant is the cost of GPU-CPU data copies! –Similarly (later): the more complex the process, the less relevant is the overhead from scalar Fortran in madevent! –And the fewer events in flight needed to fill the GPU... •In this talk I mainly give performance numbers for complex processes like gg→tҧ tgg or gg→tҧ tggg e+e-→+- 𝐠𝐠→𝐭ҧ 𝐭𝐠𝐠 A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 48 / 26 CUDACPP vs PFs - GPU ME throughputs (standalone application) •The performances of the SYCL and Kokkos implementations of MG5aMC seem comparable to direct CUDA –Further comparisons are in progress, performance scales differently with more jets for different backends (next slide) •SYCL and Kokkos run out of the box also on AMD and Intel GPUs –They also run out of the box on CPUs (performance under investigation) Xe-HP is a software development vehicle for functional testing only - currently used at Argonne and other customer sites to prepare their code for future Intel data centre GPUs XE-HPC is an early implementation of the Aurora GPU INTEL NVIDIA AMD (gg_ttgg) 16k Fixed GPU-grid size (throughput plateau) Variable GPU-grid size (throughput scan) A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 49 / 26 CUDA vs SYCL on NVidia A100 •SYCL and CUDA implementations have ~similar performances but –SYCL seems better for less complex processes –CUDA seems better for more complex processes •These are very recent results, which are still being digested (WIP!) –It will be very interesting to understand more in detail what goes on We plan to also compare more systematically the CUDACPP and SYCL performances on CPUs (vectorization, multi-core), but it will take time and optimization tweaks... WIP! PRELIMINARY! N. Nichols, T. Childers (SYCL) J. Teig (tests/plots) 𝐠𝐠→𝐭ҧ 𝐭 𝐠𝐠→𝐭ҧ 𝐭𝐠 𝐠𝐠→𝐭ҧ 𝐭𝐠𝐠 𝐠𝐠→𝐭ҧ 𝐭𝐠𝐠𝐠 CUDA < SYCL CUDA < SYCL CUDA > SYCL CUDA > SYCL A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 50 / 26 Benchmarking –Madgraph and the HEP-SCORE project •HEPscore: the new HEP benchmark for compute resources, replacing HepSpec06 –Based on reproducible HEP workloads (GEN, SIM, DIGI, REC...) within docker containers –The first version HEPscore23 should become production in April 2023 for (x86 and ARM) CPUs •The aim is to benchmark a fully loaded server: all CPU cores, and eventually all associated GPUs –(and ideally measure how well an application is doing compared to the theoretical power of the server...) –fill all CPU cores by a combination of application multi-threading and/or several identical copies/processes •A first container based on our Madgraph-on-GPU has been prepared –Very useful because it gives the same physics results on CPU and GPU: may compare them to each other! –And eventually may be used to evaluate heterogeneous processing on CPU+GPU... •The plots on the next slides are based on this HEPscore container: several identical copies/processes –(A multi-threaded CUDACPP version exists but not optimized yet –SYCL and Kokkos also provide MT) A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 51 / 26 gg→𝒕ҧ 𝒕gg (float) gg→𝒕ҧ 𝒕gg (float) ME throughput in C++ for gg→tҧ tgg (on all the cores of a CPU) •Previous tables for SIMD speedups on C++ were for a single CPU core •Large SIMD speedups are also confirmed when all CPU cores are used –AVX512/zmm speedup of x16 over no-SIMD for a single core slightly decreases to ~x12 on a full node (clock slowdown?) –Overall speedup on 32 physical cores (over no-SIMD on 1 core) is around 280 (maximum would be 16x32=512) –Aggregate MEs throughput from many identical processes using the standalone application (HEP-workload Docker container) A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 52 / 26 Some ideas for heterogeneous processing To further reduce the relative overhead of the scalar Fortran MadEvent - parallelize it on many CPU cores? •Blue curve: one single CPU process using the GPU –For gg→𝑡ҧ 𝑡gg, you need at least ~16k events to reach the throughput plateau •Yellow, Green, Red curves: 2, 4, 8 CPU processes using the GPU at the same time –Fewer events in each GPU grid are needed to reach the plateau if several CPU processes use the GPU –The total Fortran RAM would remain the same, but the CPU time in the Fortran overhead would be reduced –(Why total throughput increases beyond the nCPU=1 plateau is not understood yet!...) Throughput variation as a function of GPU grid size (#blocks * #threads) This is the number of events processed in parallel in one cycle Nvidia V100 GPU Silver 4216 4-core CPU A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 53 / 26 Lockstep beyond event-level parallelism •Efficient data parallelism (lockstep processing) requires the same function computed for different data –This is true in MG5AMC at the event level (different events i.e. different phase space points) –But it is also true at the sub-event level (different helicities within the same event) •We are evaluating the move to a different data parallelism strategy on GPUs –Currently: one event (sum over all helicities) per GPU thread –In the future: one helicity of one event per GPU thread? •Advantages: –You can fill the GPU with much fewer “events in flight” – more balanced sampling/integration in MadEvent –This is a prerequisite for moving the color matrix to externally-launched cuBLAS and tensor cores –This is also a prerequisite if we want to evaluate much smaller kernels •From all Feynman diagrams in one kernel to one Feynman diagram per kernel? •Which might decrease register pressure and increase kernel occupancy, but would require more global memory access A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 54 / 26 Beyond Madgraph •ANY matrix element event generator is a perfect fit for vectorization and GPUs! –Opportunities for 100% lockstep processing at the event-level and event the sub-event level –Had some discussions with the Sherpa team about vectorization... •Beyond matrix element event generators: are parton showers a good fit for vectorization and GPUs? –More stochastic branching than in ME event generators –But the calculations before/after each PS splitting are the same? Would just need some basketization... –Some discussions a few months ago with the Pythia team on possible collaborations A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 55 / 26 Acknowledgements •We gratefully acknowledge the computing resources provided by the Joint Laboratory for System Evaluation at Argonne National Laboratory, the Jülich Supercomputing Centre at Forschungszentrum Jülich, the Super Computing Applications and Innovation department at CINECA, and the EuroHPC Joint Undertaking at CSC’s Kajaani data centre •Many thanks from me and the whole team to: –Our CERN IT colleagues (especially Ricardo Rocha and Ulrich Schwickerath) for their help with GPU nodes! –Sebastien Ponce, Hadrien Grasland, Steve Lantz, Marco Clemencic for their suggestions on vectorization –Igor Vorobtsov and Klaus-Dieter Oertel for useful discussions on vectorization on Intel CPUs –The organizers and our mentors at the Sheffield 2020 and Lugano 2022 GPU hackathons –Domenico Giordano and the HEPiX benchmarking WG –The organizers of this Compute Accelerator Forum –The HSF event generator WG –The original authors of Madgraph5_aMC@NLO •Many thanks from me personally to the whole madgraph4gpu team for the great collaboration! ☺ A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 56 / 26 Functional tests –CI, standalone, madevent •Functional tests in the CI –Starting from some hardcoded momenta, reproduce some hardcoded matrix elements (within precision) –A few different processes, architectures, all in double-float-mixed precision •Functional tests for standalone tests –Always use the same random seeds, visually check that results do not change in logs •Functional tests for madevent –Compute cross sections and generate LHE event files in Fortran –Then do the same in CUDA/C++ and compare results (within precision) –Note: for floats, the tolerance for allowed differences is now very generous... A. Valassi –Lessons learnt from vectorizing Madgraph5_aMC@NLO ECFA Workshop, Bruxelles, 22 June 2023 57 / 26 Build challenges and issues •The build times increase enormously as you add more final state gluons –Many more diagrams... •...but also the result of some whole-program link time optimization? (we disable RDC in CUDA) –This may improve if/when we go to smaller kernels –We use ccache to ease the pain •Conversely, ninja would probably not help as we have few very large compilation units, not many small ones •We currently build separately for our five SIMD modes –Eventually we might move to a “fat binary” approach with dynamic choice of the best implementation