WIP: draft of multi-threaded `modified_equation`

[ranocha]

I have implemented a rough first sketch of a multi-threaded version of modified_equation.

@ketch This may be interesting for us but I would like to get your feedback on a possible interface first.

Currently, I extend the interface modified_equation(series) to modified_equation(series_integrator, thread = Threads.nthreads() > 1). If thread == true, it will use a multi-threaded version - which is the case by default if Julia is started with multiple threads. This allows us to check both code paths in CI without having to use multiple threads and allows users to restrict this function to the serial version even if using multiple threads.

However, I haven't extended this optional argument to all versions - something like modified_equation(A, b, c, order) will still create the B-series of the Runge-Kutta method specified by A, b, c first up to order and then use call ``modified_equation(series_integrator, thread = Threads.nthreads() > 1)`.

Preliminary results

On my system, I currently get

julia> using BSeries, BenchmarkTools

julia> series = bseries(AverageVectorFieldMethod(), 11);

julia> modified_equation(series);

julia> @benchmark modified_equation($series)
BenchmarkTools.Trial: 3 samples with 1 evaluation.
 Range (min … max):  1.744 s …  1.762 s  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.754 s             ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.754 s ± 9.139 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █                              █                       █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  1.74 s        Histogram: frequency by time        1.76 s <

 Memory estimate: 601.44 KiB, allocs estimate: 46.

julia> Threads.nthreads()
1

on current main and

julia> using BSeries, BenchmarkTools

julia> series = bseries(AverageVectorFieldMethod(), 11);

julia> modified_equation(series);

julia> @benchmark modified_equation($series)

BenchmarkTools.Trial: 11 samples with 1 evaluation.
 Range (min … max):  478.578 ms … 492.619 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     488.493 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   486.350 ms ±   5.482 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █ █ █      █                             █ █   █  █ █   █   █  
  █▁█▁█▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁█▁▁▁█▁▁█▁█▁▁▁█▁▁▁█ ▁
  479 ms           Histogram: frequency by time          493 ms <

 Memory estimate: 618.31 KiB, allocs estimate: 229.

julia> 

julia> Threads.nthreads()
4

with this PR.

TODO

• [ ] Clean up internal code (see TODO notes)
• [ ] Extend approach to modifying_integrator
• [ ] Add tests
• [ ] Add benchmarks?
• [ ] Add multi-threaded CI job?
• [ ] Fix multi-threaded buffers in RootedTrees.jl, cf. https://julialang.org/blog/2023/07/PSA-dont-use-threadid/