Performance Engineering master plan

[mmesiti]

List of aspect to tackle, in order of increasing estimated cost/benefit ratio:

1. why workspace does assignment take so much time?
2. why is the solver taking so much time? Likely because of processing related to callbacks [TO CONFIRM]
   1. try different integrators. Result: the % of time spent outside of getDeriv! does not actually change.
   2. parallelization of solver update with MPI? <- this is likely a bad idea, solver updates are likely cheap and communication between nodes would make them even slower?
   3. make sure threading is used in the status update? update is good, callbacks might be at fault here
   4. optimize/tweak callbacks. Asynchronous callbacks?
3. write a function that returns the is, it and iu ranges in a more balanced way with a bias towards splitting iu and it ranges first, plug that into the current code and see what happens (This is not optional, the current partitioning logic is flawed, but if only the s-dimension is split it is not a problem).
4. keep experimenting with the custom data structure (see Optimized data structure for X.{a,b,c} #2 ) and see if there are improvements (it makes MPI parallelization easier, but the custom {set,get}index[!] functions might (but should not) cause performance issues.
5. keep experimenting with code generation

[mmesiti]

Regarding the solver update: I this is not done in parallel with Threads. When I change the number of threads, for example, I see (using BS3, N=25 and Nlen=14):

nthreads per rank	time for solver	time for deriv!
76	12.1s	8.6s
38	18.2s	14.8s
19	30.4s	27.0s
10	51.4s	48.1s
5	97.2s	93.9s
2	232s	229s
1	459s	455s

And very interestingly, the time for the solver is always 3-4 second more than the time for deriv! for all number of threads, meaning that

1. either it is tragically memory bound, or the problem is to the point that 1-2 threads saturate the performance
2. or it's not multithreaded
   My suspicion is that lack of multithreading might be the issue. Every solver for ODE should be doing many times operations of this kind:

state .+= deriv .* dt

that should benefit from multithreading to some extent. If I understand the situation correctly, we are using ArrayPartitions for state and deriv, and although RecursiveArrayTools provides all the methods required by Base on AbstractVector, I do not see any sign of multithreading being used.

Trying this code:

julia> A = rand(50,50,50);

julia> B = rand(50,50,50);

julia> @btime B .+=  A .* 4.5;
  76.862 μs (4 allocations: 128 bytes)
  
julia> @btime @tturbo B .+=  A .* 4.5;
  14.159 μs (4 allocations: 128 bytes) 

Meaning that for a 50^3 array we should see at least a factor of 5 when using threading (in the limit for high thread numbers).
So I guess we might try to provide the methods for StateType that are needed by solve and make sure threading is used there, instead of wrapping it in an ArrayPartition?

[mmesiti]

Continuing research on this idea that the inner logic of the solver is not multithreaded, I tried to find with methods are used by the solver internally, by profiling the code again. It is my understanding that one of the important parts could be the call to the interpolant, here (at least, this is in the output of the profiler). Do I read the code correctly? It does not seem multithreaded - am I missing something?

[NilsNiggemann]

yes, that might explain things, (although it should be multithreaded)!
Between steps, the solver calls a "SavingCallback" (and a probably negligible callback for output printing). This callback computes observables from the full ODE state that are to be saved during the flow. So possibly, the problem is one of the following:

1. Since the time steps at which the observables are computed do not correspond to ODE steps, the library will use hermite interpolation of the ODE solution. Perhaps constructing this interpolation is more expensive than I thought? -> This might explain why the overhead does not decrease with the number of threads.

2. The relevant part of the callback itself. Although this is less demanding than the ODE rhs, this could in principle also become somewhat expensive. Heres the function that might be the culprit:
   getChi
   This is multithreaded though, so in principle it should scale relatively well.

[NilsNiggemann]

If the problem is the construction of the Hermite interpolation, we might have a tough problem. For testing, we should be able to decrease the overhead with the option SolveFRG(Par, ObsSaveat = [...]) where we provide a custom array of steps in the ODE integration, where the observables should be saved.

Another option for testing would be to remove the saving callback and instead move the "savefunc" to be called in the functionCallingCallback instead.
This would ensure that the function is only called once at the end of each step so there would be no hermite interpolation. Of course this would be a problem for practical application of the code where it is desirable to have smooth values of the observables, so it would only be for testing.
Let me know if this would help, I could make a quick PR!

[mmesiti]

Sorry for taking so long to reply. These are the times I measure with TimerOutputs.jl, when using the BS3 solver and N=25, on 2 nodes:

────────────────────────────────────────────────────────────────────────────────────────────────────────
                                                    Time                                Allocations
                               ───────────────────────────────────────────────   ────────────────────────
       Tot / % measured:                        12.4s /  97.3%                        441MiB /  80.8%

 Section               ncalls     time    %tot     avg     min     max     std     alloc    %tot      avg
 ────────────────────────────────────────────────────────────────────────────────────────────────────────
 total solver               1    12.0s  100.0%   12.0s   12.0s   12.0s    NaNh    357MiB  100.0%   357MiB
   getDeriv!               40    8.64s   71.7%   216ms   204ms   244ms  8.63ms   42.1MiB   11.8%  1.05MiB
     getXBubble!           40    8.36s   69.4%   209ms   197ms   237ms  8.59ms   29.4MiB    8.2%   752KiB
       partition           40    7.91s   65.6%   198ms   187ms   225ms  7.05ms   27.6MiB    7.8%   708KiB
       communication       40    452ms    3.8%  11.3ms  9.41ms  44.7ms  6.05ms   1.26MiB    0.4%  32.3KiB
       get_ranges          40   2.63ms    0.0%  65.7μs  58.1μs  74.3μs  3.10μs    479KiB    0.1%  12.0KiB
     workspace             40    114ms    0.9%  2.85ms  2.59ms  3.12ms   143μs   2.95KiB    0.0%    75.6B
       setZero - X         40   79.5ms    0.7%  1.99ms  1.73ms  2.18ms  95.4μs      672B    0.0%    16.8B
       setZero - Deriv     40   33.7ms    0.3%   844μs   730μs   947μs  71.8μs     0.00B    0.0%    0.00B
       basic constr...     40   30.2μs    0.0%   755ns   594ns  1.03μs   131ns     0.00B    0.0%    0.00B
     symmetrizeBubble!     40   62.9ms    0.5%  1.57ms  1.03ms  3.19ms   607μs   3.69MiB    1.0%  94.4KiB
     addToVertexFro...     40   53.3ms    0.4%  1.33ms  1.22ms  1.76ms   116μs   1.86MiB    0.5%  47.6KiB
     get_Self_Energy!      40   22.5ms    0.2%   563μs   356μs  1.01ms   198μs   4.89MiB    1.4%   125KiB
     symmetrizeVertex!     40   17.1ms    0.1%   428μs   273μs   503μs  67.7μs   2.32MiB    0.7%  59.4KiB
     getDFint!             40    120μs    0.0%  3.00μs  2.48μs  3.44μs   274ns     0.00B    0.0%    0.00B
     setup                 40   61.2μs    0.0%  1.53μs   967ns  2.32μs   309ns     0.00B    0.0%    0.00B
   get_observables        178    481ms    4.0%  2.70ms  2.24ms  3.85ms   408μs   9.34MiB    2.6%  53.7KiB
   bareOutput               9    322ms    2.7%  35.8ms  2.69μs   118ms  53.8ms   6.00MiB    1.7%   682KiB
────────────────────────────────────────────────────────────────────────────────────────────────────────

Which seems to suggest that the time spent in getObservables and in bareOutput are just 4.0% and 2.7%, and getDeriv! is around 72%. This means that the time spent in the "internal logic" of the solver is ~20%. If I understand well enough what you are saying, is that in that remaining 20% most of the time might be spent to interpolate the results to compute the observables to store. What I was able to see (today) from the PProf.jl seems to confirm that, but I have tried to comment out the callbacks from the solve invocation and I was immediately confronted with errors coming from the HDF5 library, so I have not been able to verify that the overhead comes from that with 100% certainty (although, also the text output of Profile seems to confirm that the non-multithreaded part that I saw is part of the callback logic)

I also tried to set ObsSaveat = [] but I got errors which at the moment I cannot explain.

To me, any task related to a saving callback should actually be embarrassingly parallel - after copying out the state of the solver to a new object (possibly doing it multiple times if this is required by the interpolation), the computation related to the callback should be able run asynchronously. Are such callbacks available in OrdinaryDiffEq?

Anyway, for the moment I am happy that we know this, set a note on this, go forward and in case we need we can get back to this.

[NilsNiggemann]

Plan:

• Run again for S = getCubic(18,[1.]); Par = Params(S,N=50,T=0.4,accuracy = 1e-9), SolveFRG(Par,method = VCABM(), CheckpointSteps = 3, CheckpointDirectory = "/scratch/... " )` (hard) (or NLen = 15, T=0.6) easier. (less memory: method = DP5())
• Disable callbacks completely --> performance regression gone?
• Enable callbacks but callback function is trivial (test whether it is hermite interpolation)
• If it is hermite interpolation: Try other arrays instead of partition arrays

[mmesiti]

Reminder: keep further discussion here: NilsNiggemann/PMFRG.jl#8