GPU (CUDA)-kernel for end-user analysis

[Moelf]

Assuming we eventually want to do end-user analysis on GPUs, the big picture looks like:

Although histograming is not the biggest bottleneck (for example, we don't know how to do GPUDirect for .root files at the moment), it is something we have to solve:

To that end, scikit-hep has: https://github.com/scikit-hep/cuda-histogram, but at the moment, it simply calls cupy kernels:
https://github.com/scikit-hep/cuda-histogram/blob/ef17af76959683e461fff4c73d28d2778cd658e9/src/cuda_histogram/hist.py#L244

this has a few problems:

• It's not very fast (cupy is not very fast, and we launch multiple kernels to get weights and sumw2)
• it doesn't fuse the kernel with analysis loop

In the long run, in order to fuse analysis loop and histogram kernel, we probably need both of them to be written in Jax (or something).

[lgray]

The predecessor of https://github.com/scikit-hep/cuda-histogram is scikit-hep/coffea#1095

Another important thing to consider in histograms is that frequency of bin-hits across warps = slow down. So higher dimensional histograms are faster (compared to boost-histogram, even cupy is 10x faster at 4 dims and a large number of bins). At lower dim, it is actually quite easy for histograms to become the bottleneck, on GPUs. Implementation matters a lot.

There is also: scikit-hep/uproot5#1426 for storage integration. Where we are most of the way to using GPUDirect to read root files (well, RNTuple, TTree is very much an anti-pattern for GPUs).

A major thing missing is corrections on GPUs.

Not all awkward GPU kernels are very well optimized yet.

Fusing kernels comes with jit compilation cost.

numba-cuda is much more difficult to learn and understand than just bare numba, so there is a user training component as well.

[Moelf]

Another important thing to consider in histograms is that frequency of bin-hits across warps = slow down. So higher dimensional histograms are faster

this can be mitigated using per-block shared memory, I have some slides on this I can share. Atomics into global array is not too bad but we can do much better in most HEP applications.

[lgray]

There's a nice nvidia paper on this from like.. 8 years ago as well. I'll find and post.

[Moelf]

even cupy is 10x faster at 4 dims and a large number of bins)

the problem of this is we're not considering a system with 1 GPU usually has more than 10 CPU cores. I currently have an implementation that's 270x faster asymptotically (after ~5*10^4 input length) and I still consider it to be too slow, because we can easily find a system with many cores these days.

There's a nice nvidia paper on this from like

Nvidia's reference histogram example never got updated to use shared memory, but they also don't consider weights and sumw2 (I suppose they don't matter much, just more data movement...).

[lgray]

A large class of people do not care about weights or sumw2, but it is getting smaller as more disciplines realize how useful carefully considered weighted data can be.