API (5) Encapsulate the ME calculation in a single class #301

This is a placeholder for one of the next big steps in the development process: an API to encapsulate the ME calculation in a single class. Let's call this a "ME engine" or "ME kernel launcher" for the time, say MEKL.

The idea is to make it easy to plug and play the ME calculation into different types of workflows, these two in particular:

1. A standalone C++/CUDA application as we have it now, specifically:

(1a) CPU only, pure C++. A sampler produces hstMomenta on the host. The MEKL produces hstMEs on the host.

(1b) GPU, in C++/CUDA. A sampler produces devMomenta on the device. The MEKL produces devMEs on the device. Someone (probably another component, not MEKL?) may also need to copy devMomenta to hstMomenta on the host and/or devMEs to hstMEs. But neither of these two extra tasks needs to necessarily happen all the time: you may decide to only compute cross sections and not need to copy momenta to host as hstMomenta (for writing to file) at all, and you may want to aggregate devMEs into avreage cross sections on the device, without copying them to hstMEs at all.

2. An external application such as the MadEvent-based Fortran application, specifically:

(2a) CPU only. A sampler like MadEvent produces momenta on the host. These may be hstMomentaFortran in Fortran matrix convention, which may need a transposition to hstMomenta in the C++ matrix convention needed by the MEKL. The MEKL produces the hstMEs: this is the same components as in case (1a) above.

(2b) GPU. A sampler like MadEvent produces momenta on the host. These may be hstMomentaFortran in Fortran matrix convention, which may need a transposition to hstMomenta in the C++ matrix convention needed by the MEKL. Then someone copies the hstMomenta to devMomenta (or maybe it is the hstMomentaFortran that are copied to devMomentaFortran and it is then the GPU that does the transposition to devMomenta - I think this is what Stefan does now). Then the MEKL computes devMEs from the devMomenta, this is exactly the same component as the MEKL above in 1b. Then the relevant outputs (devMEs individually, or some aggregates) must be passed back to the host and back to the Fortran framework.

In all of the above, the common core component (which is what spends most of the compute time for complex processes) is the MEKL computation of hstMEs from hstMomenta (CPU) or devMEs from devMomenta (GPU). This is more or less the current sigmaKin function, with a few things to change

• before the klas5X PR Remove neppV and fptype_sv from API and decouple neppV from neppM (epochX) #300, hstMomenta was an fptype_sv* and devMomenta was an fptype*; after that PR, both will use a simple fptype* array
• however, the important thing for vectorization is that the fptype* must be reinterpreted as fptype_sv*, which only needs a minimal constraint of alignment - presently in klas5X I am using a 64byte (512bit) alignment, which is enough for all types of vectorizations up to AVX512; note that also hstMEs must be similarly aligned as some operations on hstMEs are also vectorized
• the actual AOSOA structure of hstMomenta should also have neppV a multiple of neppM to allow reinterpret casts; the best performance is in any case found if neppM=neppV as it was initially
• the current klas5X code has a variety of implementations, which may also work with non-aligned data or with neppM<neppV (eg with an AOS neppM=1); for performance reasons and to keep the code clean, my suggestion is to only support very few cases and implement sanity checks to check that the constraints are met; I suggest in particular to always assume the data are aligned, and to only support one other non default readout, eg neppM=1
• note that the bulk of the calculation is vectorized independently of the way the hstMomenta are defined (aligned or not, neppM=neppV or not), because the momenta are initially transformed and passed by vale as fptype_sv anyway to the calculations in the XXX functions, an dthen remain fptype_sv throughout
• for CUDA, the reading of the AOSOA is now done in the XXX functions; eventually it may be better to do them outside (in calculate_wavefunction, eventually in the MEKL), this is PR (DO-NOTHING) reverted tests for further separating data access from calculations in xxx #189
• given this long preamble, one idea is to have the MEKL class itself prepare the appropriate input and output buffers, or in a second step accept to be constructed from well defined input and output buffers (which eg may be not aligned or have any neppM value); one important point is that any sanity checks eg on alignment would be done only once, rather than at each computation of the computeMEs() function, which would slow it down

One important point on this MEKL class is that one could have several them in parallel. For multithreading in C++, in particular, I would find it better to have N such classes (with completely separate data buffers on different NUMA domains) for N different threads, rather than rely as we do now on OMP within each computeMEs call to split up the calculation on several threads, which is observed to be very inefficient (issue #219

Voila, this is a long brain dump... to be discussed with the rest of the team!

I am about to send a very draft prototype of a WIP PR for this.