AOS/SOA for input particle 4-momenta (and random numbers)

[roiser]

Done, including AOSOA (AV)

eemumu_AV/master

Steered by typedefs. We need to go to the end with this idea.

Should go back and time AOS… actually AOSOA5.0E8, AOS4.8E8

[valassi]

Actually I have timed this yesterday, and AOSOA does seem to pay off, even if not by much.
See https://docs.google.com/document/d/1g2xwJ2FsSlxHvSUdPZjCyFW7zhsblMQ4g8UHlrkWyVw/edit#

Description	CUDA tput
BASELINE	5.00E8/s
SOA	4.65E8/s speed -7%
AOS	4.85E8/s speed -3%

I would say that we should keep some AOSOA-like structure as baseline

To do in any case

• understand better what causes the performance improvement (memory access via coalescing? SIMD vectorized computations??)
• massive clean up of the code, should use the same array dimensions everywehere: goal should be AOSOA[ngpublocks][nparticles][np4=4(E,px,py,pz)][ngputhreadsinblock]... also the random number arrays should use ngputhreads as last dimnesions (now it is hardoded as 32, the number if threads in a warp)

[valassi]

Some food for thought from Vincenzo, found by chance from google "godbolt cuda soa"
https://indico.cern.ch/event/851670/contributions/3585184/

[valassi]

And another interesting talk at CERN from the same google for godbolt
https://indico.cern.ch/event/932905/contributions/3920347/

[valassi]

I would like to do some changes to the AOSOA, but first I's like to understand a bit better if/why the layour makes a difference. I can think of two things, memory coalescing, and instruction vectorization. So I am doing some research to understand which metrics are relevant, eg in the profiler.

About memory coalescing, google "nvidia nsight coalesce" brought me here: https://developer.nvidia.com/blog/using-nsight-compute-to-inspect-your-kernels. This is avery interesting article pointing to two metrics, and showing how to focus on them in the tools. I added them here 0429183

I have then profiled the BASELINE ASA against AOS and SOA. At face value I confiorm that SOA is similar, just ~2% slower than ASA, while AOS is quite a bit slower, around 7-10%. The profiles are very interesting.

First, looking at the metrics in the blog above, indeed the number of requests for memory in sgnakin is the same in ASA and AOS, but the number of sectors (transactions) is a factor 4 higher with AOS than with ASA. This is a clear indication that the memory (the allmomenta memory) is not coalesced in AOS and is better in ASA. The distinction between ASA and SOA is much more subtle and unclear, the numbers of requests and sectors is lower in SOA than ASA, with a ratiob etween them that remains comparable. Note also that the ratio is around 6, while the optimale should be around 4?

This is ASA

This is SOA

This is AOS

Second, SOA has a higher number of registers that ASA (182 against 152) and this may be the reason for a penalty elsehwere, eg in the number of active warps. Maybe this explains the slightly lower throughput of SOA, maybe not. TThis is SOA compared to ASA

[valassi]

I have made many studies related to this and I will dump a few results an ddecisions here.

First, I cleanly spearated in the code the ASA structure for random (based on neppR) and for momenta (based on neppM). In issue #23 I studied the impact of doing a compile time vs a runtime choice for ASA parameters. The difference is small but visible, at the levl of 5% throughput for neppM, and visble in other profiler parameters.

[valassi]

Next, I studied whether AOS is any different from ASA with neppM=1. They are actually almost identical. There are differences below 1% in some metrics and I am not even sure why.

This is an overview

This is the warp state

This is the instruction mix

The #registers are essentially the same.

Conclusion: I will drop AOS from the code to make it simpler. This can always be recovered by setting neppM=1 in the ASA option (and indeed it is interesting for some studies, see later).

[valassi]

Then, I compare AOSOA (in my default with neppM=32) to SOA (which is essentially AOSOA with a much larger neppM=16384*32). The latter is worse in all relevant metrics an ddoes give a lower throughput. I am not sure why, but the issues seem to come from a much larger number of registers (operations not vectorized??). Anyway, even from first principles this is not really a sound choice.

Overview

Memory

Compute workload

Scheduler and stalls

Insruction mix

Conclusion: I will finally drop also AOS. I will then concentrate only on AOSOA. This will allow a much cleaner code. The relevant parameter in any case is neppM (smaller as 1 means AOS, larger as ndim means SOA) and various options can still be studied.

[valassi]

Coming now to the optimal nepp size in AOSOA. So far I am using a default of 32, which I included as this is the number of threads in a warp. I am instead getting the idea that 4 is a better size for our doubles (8 bytes), because the cache lines are 32 bytes and fit 4 doubles. (Probably 8 would be better for floats).

To start with, at face value an ASA32 vs ASA04 (where npeeR=neppM=32 or 4). The throughput is marginally better. Especially, the FP64 seems a bit better utilized, and in the end that's the only thing we do: compute a lot of FP64, so we should do that as fast as possible.

Note that the memory difference comes from a large increase in data pipe LSU wavefront, but I am not sure that's any good in itself. I would say 4 is better than 32 because it takes overall a bit less time, and FP64 is better (both only at the level of 1%).

I also compared ASA04 to ASA08, but the latter is not better, and maybe slightly worse, for doubles.

Instead, ASA01, ie AOS, is very interesting. This is clearly worse than ASA32 (or ASA04). The big difference here is that the number of requests is the same, but the number of sectors is a factor 4 higher: this is an indication of non-coalesced memory access.

For the record I also tried ASA02, but clearly this needs twice as many transactions as ASA04.

Note that indeed for FLOAT the minimum top get fully coalesced access is ASA08 (8 floats of 4 bytes are a cache line of 32 bytes). Using ASA04 there results in twice the number of transactions as one could use.

Conclusion: in other words, we need at least ASA04 (double) or ASA08 (float) to have coalesced memory access. We do not seem to get anything better by using higher nepp. I will set the new defaults.

[valassi]

Next point is to try and understand the number of requests. I would imagine that to improve on memory usage we can do two things:

• First, for a given number of requests, try to serve them with the least number of transactions. This is about reducing the transaction/request ratio, and is what I did above. As long as we use at least nepp=4 for doubles and nepp=4 for floats, then access is coalesced and we minimise the transactins neede for a given request.
• The next point is try to understan the number of requests themselves.

To study this I used the two metrics described here https://developer.nvidia.com/blog/using-nsight-compute-to-inspect-your-kernels/. I actually realised that they are available also in command line, eg as

/usr/local/cuda-11.0/bin/ncu --metrics l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum,l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum ./gcheck.exe -p 16384 32 1
...
  gProc::sigmaKin(double const*, double*), 2020-Aug-17 17:09:35, Context 1, Stream 7
    Section: Command line profiler metrics
    ---------------------------------------------------------------------- --------------- ------------------------------
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum                                request                      1,527,808
    l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum                                  sector                     10,154,052
    ---------------------------------------------------------------------- --------------- ------------------------------

The question is why 1527808 requests? And why 10154052 transactions by the way?

Using my current default nepp=4, I ran this for several configurations

-p blocks threads iterations	requests	sectors (transactions)
-p 16,384 32 1	1,527,808	10,154,052
-p 2,048 32 1	352,256	2,558,406
-p 256 32 1	70,144	528,896
-p 32 32 1	8,768	66,112
-p 2 32 1	548	4,132
-p 1 64 1	548	4,132
-p 1 32 1	274	2,066
-p 1 4 1	274	274

This is a very useful link: https://stackoverflow.com/questions/60535867. I actually found it by looking fo rthe names of the old metrics (https://docs.nvidia.com/nsight-compute/2019.5/NsightComputeCli/index.html#nvprof-metric-comparison) because the new ones are still not much documented.

A few comments on teh above:

• The first point is that each request is a request at the warp level, which is 32 threads. Whether I use 4 threads or 32, the number of requests is the same, 274.
• The second point is that a transaction is a 32-byte cache line. So, one coalesced request for doubles (8 bytes) at the warp level needs 8 pages of 4 doubles. However, other types of requests may need a differnt number of transactions. Here it seems that 274 is 256+18. These requests for a full warp need 2066 transactions, which is 256x8(=2048)+18. So, maybe this is 256 doubles, and 18 something else.
• As soon as I request more warps (either with more threads than 32 per block, or with more blocks), the requests increase. Both for one 64-thread block, or two 32-thread blocks, I need two warps, and twice the number of requests. The number of transactions also doubles.
• This linearity goes on until 256 blocks: "256x32" needs 256 as many requests and transactions as "1x32", 70144=256x274 and 528896=256x2066. For larger numbers of blocks I cannot really understand the pattern anymore, the number of requests decreases with respect to what I would expect. To be understood further..

[valassi]

It turns out that most of the 274 requests were due to the way I had implemented helicity filtering. I have improved helicity filtering (issue #24) and this now reduces drastically the number of requests. For the default nepp=4,

./profile.sh -nogui -p 1 4 1
  gProc::sigmaKin(double const*, double*), 2020-Aug-18 09:25:37, Context 1, Stream 7
    Section: Command line profiler metrics
    ---------------------------------------------------------------------- --------------- ------------------------------
    l1tex__t_requests_pipe_lsu_mem_global_op_ld.sum                                request                             16
    l1tex__t_sectors_pipe_lsu_mem_global_op_ld.sum                                  sector                             16
    ---------------------------------------------------------------------- --------------- ------------------------------

I redo my previous table with the new code

-p blocks threads iterations	requests	sectors (transactions)
-p 1 4 1	16	16
-p 1 32 1	16	128
-p 1 64 1	32	256
-p 2 32 1	32	256
-p 256 32 1	4096	32768
-p 2048 32 1	32768	262134 -- 262141
-p 16384 32 1	262144	2097044 -- 2097063
-p 2048 256 1	262144	2097152

The number of requests follows a perfect linearity: with 16384/32 I need 16384 as many requests (262144=16384*16) as with 1/32. In other words, each warp of 32 threads issues 16 requests, and it's not surprising as it needs 16 doubles (the 4-momenta of 4 particles).

The number of transactions also seems to follow the same linearity. For high number of blocks the number of transactions fluctuates (maybe the profiling tool is not very precise on this point??), but essentially with 16384/32 I also need 16384 as many transactions as with 1/32 (as 16384*128=2097152... this is actually what I get with 2048/256, but the throughput is smaller than with 16384/32).

In particular

• With 1/4 ie only 4 threads each request needs only one transaction: each request is for one double, ie 8 bytes, but there are only 4 threads in a warp so each warp gets excatly one page of 4 doubles in a 32-byte transaction.
• With 1/32 ie 32 threads, each request needs eight transactions (I need eight 32-byte pages to serve the 32 doubles I need for one request)
• NB one "request" is for instance the request for pz of particle 2, if I understand correctly

This looks much better.

[valassi]

I have now repeated some of the previous profile tests on ASA with nepp=4 (default now) vs 32 (old default) vs 1 (ie AOS).

I would say that ASA32 and ASA04 are indistinguishable, while ASA01 ie AOS is different. Surprisingly, it seems to have now a bit higher throughput?

Anyway the memory is clearly less optimized: it is much more used because the data is retrieved in too many transactions. So there is a 75% L1 hit rate with AOS/ASA01, while it is 0% with ASA04, simply because each page is retrieved 4 times, so 3 times out of 4 it is already in cache... But it is much better to just retrieve it only once

I would CONCLUDE this issue on AOS/SOA structures FOR MOMENTA

• It does not seem to be the primary issue: our main issue is that we must optimise the internals of sigmakin, not the retrieval of the momenta into sigmakin
• That said, it seems reasonable to keep nepp=4 (for doubles, or nepp=8 for floats) as the baseline. This reduces the number of transactions, avoids relying on the L1 cache by simply retrieving all relevant data in one go.
• The fact that there seems to be a small improvement in throughput with AOS is surprising, could be better studied, but I would just ignore it. Also, it may work here for eemumu where there are exactly four particles, but it may give issues elsewhere. Actually one particle 4-momenta is one cache line with doubles (4 doubles), but it becomes half a cache line with floats (4 floats is 16 bytes), and with odd number of particles it may give issues.
• Having the code already set up (with a variable nepp) as AOSOA may also be useful to exploit SIMD on the CPU.
• The structure of memmory layouts for intermediate results of the sigmakin computation (the w[5][6]) remains instead something that we may need to study better. THis is discussed in issue Memory layout (shared/global/local, AOSOA...) for intermediate wavefunctions in ME calculations #7. I would close this one now...

[valassi]

Marking as closed in the two attached projects.

[valassi]

Hi @roiser @oliviermattelaer @hageboeck I did a few additional tests of AOSOA in cuda while writing the paper (to avoid making stupid statements). THis updates the observations above. This is all in PR #209, just merged.

Two main observations:

• I confirm that in CUDA eemumu the AOSOA for momenta only brings a very marginal improvement with respect to AOS. I do see that using an AOS moves around four times more memory (the number of sectors goes up by x4), but the ME throughout only decreases by maybe 5%. And this is eemumu, where memory access to momenta is much more important than in ggttgg...
• Second point, the number of requests that is now needed per event (or better, per event page) is 40: we need the full 4-momenta of two particles, the pz momenta of two particles (thanks to the simplified ixx/oxx functions), all this times four good helicities. In my previous posts here I was mentioning 16, but this was for a single helicity, and before the simplification in ixx/oxx.

So, all looks understood here...

This is the current default
214b8c2

On itscrd70.cern.ch [CPU: Intel(R) Xeon(R) Silver 4216 CPU] [GPU: NVIDIA Tesla V100S-PCIE-32GB]:
=========================================================================
Process                     = EPOCH1_EEMUMU_CUDA [nvcc 11.0.221]
FP precision                = DOUBLE (NaN/abnormal=0, zero=0)
RanNumb memory layout       = AOSOA[8] [HARDCODED FOR REPRODUCIBILITY]
Momenta memory layout       = AOSOA[4]
EvtsPerSec[MatrixElems] (3) = ( 7.227665e+08                 )  sec^-1
EvtsPerSec[MECalcOnly] (3a) = ( 1.352993e+09                 )  sec^-1
MeanMatrixElemValue         = ( 1.371706e-02 +- 3.270315e-06 )  GeV^0
TOTAL       :     1.119368 sec
     3,329,955,135      cycles                    #    2.649 GHz
     4,782,880,833      instructions              #    1.44  insn per cycle
       1.414098430 seconds time elapsed
==PROF== Profiling "sigmaKin": launch__registers_per_thread 120
==PROF== Profiling "sigmaKin": sm__sass_average_branch_targets_threads_uniform.pct 100%
==PROF== Profiling "sigmaKin": requests 40 sectors 80 [-p 1 8 1]
==PROF== Profiling "sigmaKin": requests 40 sectors 320 [-p 1 32 1]
==PROF== Profiling "sigmaKin": requests 655,360 sectors 5,243,035 [-p 2048 256 1]
=========================================================================

This is with an AOS for momenta (neppM=1... I still keep random numbers as AOSOA[8] to get the same physics results)
380f06c
You can see that the number of sectors (transaction roundtrips) has gone up by a factor 4 because memory is not coaleseced. Still, thoughput only goes down from 7.20E8 to 7.00E8 (and there are huge fluctuations on this machine).

On itscrd70.cern.ch [CPU: Intel(R) Xeon(R) Silver 4216 CPU] [GPU: NVIDIA Tesla V100S-PCIE-32GB]:
=========================================================================
Process                     = EPOCH1_EEMUMU_CUDA [nvcc 11.0.221]
FP precision                = DOUBLE (NaN/abnormal=0, zero=0)
RanNumb memory layout       = AOSOA[8] [HARDCODED FOR REPRODUCIBILITY]
Momenta memory layout       = AOSOA[1] == AOS
EvtsPerSec[MatrixElems] (3) = ( 6.997163e+08                 )  sec^-1
EvtsPerSec[MECalcOnly] (3a) = ( 1.272792e+09                 )  sec^-1
MeanMatrixElemValue         = ( 1.371706e-02 +- 3.270315e-06 )  GeV^0
TOTAL       :     1.208810 sec
     3,334,791,948      cycles                    #    2.652 GHz
     4,801,124,639      instructions              #    1.44  insn per cycle
       1.502558668 seconds time elapsed
==PROF== Profiling "sigmaKin": launch__registers_per_thread 120
==PROF== Profiling "sigmaKin": sm__sass_average_branch_targets_threads_uniform.pct 100%
==PROF== Profiling "sigmaKin": requests 40 sectors 320 [-p 1 8 1]
==PROF== Profiling "sigmaKin": requests 40 sectors 1,280 [-p 1 32 1]
==PROF== Profiling "sigmaKin": requests 655,360 sectors 20,971,817 [-p 2048 256 1]
=========================================================================

[valassi]

PS Uh! I am silly.

I forgot to say, I added a NEW metric, where I do NOT include the device to host copy.

In the numbers above, the difference between AOSOA and AOS ia bit bigger: 1.35E9 to 1.27E9. But still, only around 8%.

I will quote in the paper that it seems below 10%...

[valassi]

PPS A few observations in PR #210 for single precision: not surprinsingly,

• with the optimal AOSOA[8], the number of transactions is half of that in double precision with the optimal AOSOA[4], for the sam enumber of request
• with the suboptimal AOS, the number of transactions is 8 times asmany as with AOSOA[8]... actually it is the same as for double precision with AOS

In other words:

• with single and double precision, if AOS is used, the number of requests is the same in the two cases (approximately 40 times the number of events, divided by 32) and the number of transactions is also the same (approximately 40 times the number of events) - a single request is issued for 32 threads, resulting in 32 transactions
• if the optimal AOSOA is used (4 for double, 8 for single), the number of transactions decreases by a factor 4, or 8, respectively - because one request for 32 threads can be handled with only 8, or 4, respectively, transactions, insteda of 32