splitting vecsize_memmax into wrap_size and nb_wrap

[oliviermattelaer]

status of the implementation

• update python interface to be in control of the new parameter
• update the vector.inc to have the new parameter in fortran
  • VALIDATION SIMD mode should work as before (but parameter nb_wrap is ignored
• update dsample to group by wrap the symmetry selection
• update the selection_grouping to be applied on the latest subset
  • symmetry parameter still integer at this stage
  • VALIDATION SIMD mode should work as before for one wrap
  • VALIDATION SIMD mode should work as before for multiple wrap but selection is done only based on last wrap
• update iproc.imirror,icong to a vector
  • imirror_vec, iconfig_vec (NB_WRAP as size)
  • iproc stay scalar
  • pass only one value within the bridge
  • VALIDATION SIMD mode should work as before for one wrap
  • VALIDATION SIMD mode should work correctly for fortran mode
  • VALIDATION SIMD mode should work correctly for CPP/CUDA mode (but fixing one of the config)
• update IB handling (do not move it to a vector since IMIRROR is already a vectory)
• update cpp/bridge to handle multiple value for the channel information (choice: channel information passed to the bridge PER EVENT (not by wrap, so we can choose the wrap size freely)
  • VALIDATION SIMD mode should work correctly for CPP/CUDA mode (but fixing one of the config)
• check the possibility to define a nb_wrap_used

Other Validation

• check that nb_wrap can be updated via the run_card
• check that nb_wrap can be updated at output stage

[oliviermattelaer]

point to check so far

• change (python signature for the function that create the PDF line)
• vector.inc will have issue with the line
• parameter vexsize_memmac need to be added to the run_card has system parameter to allow the vector.inc to be written correctly
• stays with a single index:
  • all_G (use vecsize_memmax)
  • couplings
  • update_as_param (only use vecid)
  • update_as_param2 (only use vecid)
    Some fortran action that needs to be done
• select_grouping sginature is now set as imirror, iproc, iconf, all_wgt, WRAP_SIZE, iwrap)
• prepare_grouping_choice signature is kept intact
• update to a vector (size can be either full or nb_wrap --better but more complex?--)
  • imirror,
  • iproc (technically not needed since force to one)
  • iconf

[oliviermattelaer]

@roiser: I have tested a scan on nb_wrap with the current fortran implementation:
(here wrap_size=32)

grid size	nb_wrap	cross-section (p p > e+ e-)	timing survey	timing refine
32	1	844.9 +- 1.978 pb	1s	12s
320	10	842.8 +- 1.973 pb	2s	10s
3200	100	844.6 +- 1.055	4s	34s
8000	250	843.2 +- 1.028 pb	5s	47s
16000	500	902.9 +- 31.41 pb	7s	69s
32,000	1000	841.9 +- 1.274	16s	0s
320,000	10,000	776.9 +- 52.13	69s	0s
3,200,000	100,000	not enough memory

So this worries me a bit, looks like at some point we have an issue like overfitting (so maybe this is also related to the fact that I pick a simple process).
In principle this has nothing to do with the wrap_size but rather with the grid size that automatically increase the number of events by iteration (for the first few iterations).

Did you try something like that on GPU with the master branch?
So I would suggest to avoid grid size > 10k for the moment.

[roiser]

Hi @oliviermattelaer here are the numbers for the GPU (A100). I have redone your numbers in the first table and then a table with more optimized ones IIUC, with vector sizes modulo 108 (which are the number of SMs on the A100).

The seconds to my mind are not conclusive but I also counted the number of refine jobs that are run which go to 0 from a certain stage on. Could this be the same indication of overfitting?

vector size	cross-section (p p > e+ e-)	timing survey	timing refine	# refine jobs
32	846.1 +- 2.143 pb	3s	45s	14
320	842.6 +- 1.857 pb	3s	18s	13
3200	845.1 +- 1.193 pb	7s	27s	11
8000	843.7 +- 1.13 pb	8s	35s	11
16000	845 +- 0.9486 pb	13s	53s	11
32000	843.1 +- 0.9455 pb	23s	53s	6
40000	842.9 +- 1.103 pb	29s	29s	0
320,000	not enough memory

vector size	cross-section (p p > e+ e-)	timing survey	timing refine	# refine jobs
3456	844.7 +- 1.226 pb	7s	27s	11
6912	843.8 +- 1.09 pb	9s	33s	11
13824	843.2 +- 1.252 pb	13s	28s	6
27648	842 +- 1.044 pb	20s	47s	6
31104	843.6 +- 1.032 pb	22s	51s	6
34560	843.5 +- 1.243 pb	24s	24s	0
38016	842.7 +- 1.123 pb	26s	26s	0
41472	841.9 +- 1.077 pb	29s	29s	0
55296	not enough memory

[oliviermattelaer]

Those are much better number that the ones I got.
This is likely a good news meaning that the issue is likely related to my wrap implementation and not to a real over-fitting.
So I will need to investigate this to understand and ( more importantly) fix it. Thanks for the check

[roiser]

I would not give too much importance to cpu seconds though, I guess what we are mostly measuring there is the fortran part.

[oliviermattelaer]

I agree for the timing, the critical point here is the fact that you have a stable cross-section which I do not have at all.

PS: I have pushed the change for having the list of channel (with name channels --with a s--). Remember that you need to update both the plugin (new_interface_wrap) and the MG5aMC interface (gpucpp_wrap). Thanks

[valassi]

Hi @oliviermattelaer @roiser thanks for all your work on this. Just a couple of points from a discussion with Stefan, since we were trying to remember what he had agreed in the past and I have been unable to find any notes about that.

First Olivier a quick point: eventually can you please change all "wrap" into "warp" in the code? Even today I got confused again what wrap is! ;-)

Then about the main issue: my main question is what cudacpp should do internally.

• My understanding is that it is Fortran that know the warp size (I would actually hardcode this to 32 also for CPU?), ie how the channelIds need to be segmented. For instance, in the channelid array you may have 32 4's, then 32 5's, then 32 6's etc. Then cudacpp receives this array.
• I am quite sure I had made the point, and we had all agreed, that there should be no hidden assumption. I mean, cudacpp should NOT assume that there are 32 x's, then 32 y's, then 32 z's. (I had also said I think that if we do some assumption we should check it with an assert, but I am not sure this is what we want).
• I think at some point I had proposed (in a discussion with @hageboeck too?) that we can simply change this code

      if( channelId == 1 ) numerators_sv += cxabs2( amp_sv[0] );
      if( channelId != 0 ) denominators_sv += cxabs2( amp_sv[0] );

into something like

      channelids_sv = CHANNEL_ACCESS::kernelAccess( pchannelIds ); // the 4 channels in the SIMD vector
      bool_sv mask_sv = ( channelids_sv == 1 );
      numerators_sv += mask_sv * cxabs2( amp_sv[0] );
      if( pchannelIds != nullptr ) denominators_sv += cxabs2( amp_sv[0] );

In other words, cudacpp does not do any assumpton at all. In practice, thanks to how the forstran has created those arrays, then the masks will always be all 0 or al 1 in cpp, but this is not assumed. Similarly, in cuda (where _sv is a scalar in a kernel), all threads in a warp will end up having exactly the same mask value and they will all be in lockstep: but again, this depends on what fortran has made of the array, and there is no assumption in the code that the threads will be in lockstep.

• Giving the above, there is no need to pass warp_size from fortran to cudacpp... and I think we ahd agreed on this, right?

Thanks... just trying to refresh my memory. Andrea

[oliviermattelaer]

Hi,

   First Olivier a quick point: eventually can you please change all "wrap" into "warp" in the code? Even today I got confused again what wrap is! ;-)

consider it as done (just need to do a git push)

      if( channelId == 1 ) numerators_sv += cxabs2( amp_sv[0] );
      if( channelId != 0 ) denominators_sv += cxabs2( amp_sv[0] );
into something like

      channelids_sv = CHANNEL_ACCESS::kernelAccess( pchannelIds ); // the 4 channels in the SIMD vector
      bool_sv mask_sv = ( channelids_sv == 1 );
      numerators_sv += mask_sv * cxabs2( amp_sv[0] );
      if( pchannelIds != nullptr ) denominators_sv += cxabs2( amp_sv[0] );

Please put comment around those lines since you pass from an easy to read code to a code that I simply do not understand (which I'm fine with obviously but I would be nice to keep readibility whenever you have those cryptic accessor.

Next question is about performance are we not shooting ourself with such complicated accessor?
But in any case, yes that was the point.

  *   Giving the above, there is no need to pass warp_size from fortran to cudacpp... and I think we ahd agreed on this, right?

This is indeed our current strategy/agreement indeed.

But clearly it has some drawback that I quote here for reference:

1. You do transfer many useless data from cpu to gpu (since each event has it's own information rather than each warp)
2. You could have set those values within the shared_memory for GPU (reducing register pressure) --which would have been a natural fit--
3. You need to implement this masking procedure (where we hope that we will never need it in practise)

While the advantages are

1. easier interface from fortran to cpp
2. possibility to compile for multiple hardware and be bullet proof for future hardware --if the warp_size is too large to ensure that MC is not bias--.
3. possibility to not trust that fortran

[roiser]

Oh, I may have completely misunderstood you then @oliviermattelaer . Are you saying you are passing in a vector of size [n events / warp size] instead of [n events]? I remember we discussed that but I thought we had agreed on n events. Either is of course fine with me.

I have for now developed everything in the branch to pass the vector in until calculdate_wavefunctions, now I need to handle the individual array entries (so e.g. the code that we discussed in a previous meeting with masking and that @valassi has spelled out above).

[oliviermattelaer]

Hi @roiser,

Are you saying you are passing in a vector of size [n events / warp size] instead of [n events]?

I do pass a vector of size [n events].
So this is has you did remember/agreed sorry if that was not clear.

[roiser]

No I do pass in a vector of n events (NOT n events / warp size), ok then I may have misunderstood you from the comment above. All good thanks, sorry my misunderstanding, reading it again its all clear

[valassi]

No I do pass in a vector of n events (NOT n events / warp size), ok then I may have misunderstood you from the comment above. All good thanks, sorry my misunderstanding, reading it again its all clear

Thanks Olivier and Stefan!

Stefan, I think Olivier means that an alternative to passing n events would be to pass n_warp (n_events/warp_size) events. But this is not what we had agreed (we agreed to pass n events).

Olivier, I see your point, indeed we might pass less memory around if we passed n_warp channels only. The disadvantage however in my opinion is that the computation becomes more complex in cudacpp as we need to introduce that concept of warp. Yes it's debatable, maybe the performance gain would be large. Anyway, I suggest we stick to our original plan, we implement the functionality (ie the possibility to have difrent channels in different events) and then we reconsider later on?

consider it as done (just need to do a git push)

Thanks! (for wrap->warp)

Please put comment around those lines since you pass from an easy to read code to a code that I simply do not understand (which I'm fine with obviously but I would be nice to keep readibility whenever you have those cryptic accessor.

Yes good point, TODO Andrea put some comments.

Next question is about performance are we not shooting ourself with such complicated accessor?

Listen, maybe yes. Again, up for review later on... but as this is done elsewhere, I would keep the same style ebverywhere, an dthen change everything if/when we decide to change it. Personally I think that the readability is ~ok (yes I should add comments). What I find less nice is that the use of templates all over the place ends up exploding build times I fear.

My suggestion on this is: if we model the accessors differently in the future, let's follow the cpp parallel model. From what I understood of the recent seminar, the model is a bit similar, there is mdspan which takes non-owning a pointer and then some metadata for defining an accessor that depends on the chosen layout. See https://en.cppreference.com/w/cpp/container/mdspan. In principle this should eventually work for GPUs too, Anyway, also about this one lets rediscuss later ok?

[valassi]

My suggestion on this is: if we model the accessors differently in the future, let's follow the cpp parallel model. From what I understood of the recent seminar, the model is a bit similar, there is mdspan which takes non-owning a pointer and then some metadata for defining an accessor that depends on the chosen layout. See https://en.cppreference.com/w/cpp/container/mdspan. In principle this should eventually work for GPUs too, Anyway, also about this one lets rediscuss later ok?

PS As the speaker had mentioned, this has close links to kokkos or is even modeled from there. Just for reference this is a tutorial I found from kokkos about mdspan https://github.com/kokkos/mdspan/wiki/A-Gentle-Introduction-to-mdspan