Feature/meson contractions

[felixerben]

Restructuring of the MContraction::Meson module. For any combination of GammaSource / GammaSink (256 in the case of the "all" option), the old module computed the propagator product one by one. This one improves on the old module in the following ways:

• If a list of gamma pairs is given to the module (i.e. the "all" option is not used), then the input is now verified and an exception thrown for illegal input - the old one happily produced correlators with "undef" gamma matrices
• Instead of looping over each individual gamma pair, there is now an outer loop over gammaSink.
  • To make this possible, the gammaList object is sorted in the setup step so that the gammaSink are in the order they are defined in Grid
  • There is also a new temporary int-vector gammaSnk that stores which gammaSink matrices are needed
  • For each gammaSink, the product adj(q2)(g5(gSnk))*(q1) is stored in a temporary object. Then in a nested loop over the relevant gammaSource, the full propagator product to be traced over is computed.
• the module also still has the granular timers in it that measure the sliceSum and mesonConnected time. I suggest we keep this to monitor performance.

NB: This new workflow does not apply to the case of pre-sinked propagators. I am not sure whether this would speed it up and the implementation would be more complicated, as the nested loop over t needs to be accounted for.

This module is tested on my laptop, CPU, comparing the output of tests/Test_hadrons_spectrum against the version currently in develop. Using h5diff, I could verify that the output agrees up to 10e-15 in my test, where I used a solver precision of 1e-8.

One thing that I cannot explain and that puzzles me: For the (Gamma5 Gamma5) meson only, the output is the complex conjugate of the output from the module in develop. I.e. the new module produces

        DATASPACE  SIMPLE { ( 8 ) / ( 8 ) }
        DATA {
        (0): {
              0.475143,
              -2.14113e-19

where the old one did

        DATASPACE  SIMPLE { ( 8 ) / ( 8 ) }
        DATA {
        (0): {
              0.475143,
              2.14113e-19

This might just be some oddity.

The timings on my laptop are not reliable enough to say with certainty whether this module executes faster. We will have to see in the production runs in tursa

[aportelli]

Hi @felixerben, thanks. The safer check of arguments is good. However the implementation relies on the underlying integer representations of Gamma::Algebra which is too fragile as it is subject to changes. Additionally I do not see why this is necessary. You could for example use a std::map<Gamma::Algebra, std::vector<Gamma::Algebra>>, which maps each source to a list of sinks, and that can be produced by parsing the input list of pairs. If such map is called ssMap, thenthe computational loop can be

for (auto &ss: ssMap)
{
  Gamma::Algebra source = ss->first;
  // compute first part of trace, the source is ss->first
  for (Gamma::Algebra &sink: ss->second)
  {
    // second part of trace
  }
}

I hope I am not overlooking something, but a structure like that should not rely on any ordering and any integer representation of Gamma::Algebra.

Finally the complex conjugate results might just be FP rounding, but do inspect more number to be sure.

I am aware we re pressed with time. Since this module is critical to many runs I would prefer the structure is cleaned as above to go in develop, but you should go ahead running from your fork if there is not time for that.

[aportelli]

Hi @felixerben thanks! Just to check: is this ready to go? As retro-compatibility been checked against the previous version?

[felixerben]

Yes this is checked against the old version and ready to be merged.

[aportelli]

Thanks!