lgarrison
commented
3 years ago I'll note that I think an even faster version would be possible if the threads were doing sequential writes and random reads, instead of random writes and sequential reads. But the sequential write version requires a fast, parallel argsort of the cell indices. I couldn't find one readily available in C, and I think the current version is fast enough that gridlink is unlikely to be a bottleneck.
Also, astropy bot appears to be afk...?
manodeep
pllim
Hey @manodeep, here's a stab at a parallel implementation of gridlink! It was motivated by a recent example that popped up in my research where the actual pair counting was very fast (due to low sample density and low rmax), but gridlink took about the same amount of time as the counting, due to the large Nmesh, and thus large number of mallocs and single-threaded particle writes. This version does one malloc for the whole grid (in detail, one malloc for each X/Y/Z/W column of particle data), calculates cell assignments in parallel, then uses that information to do parallel writes of the particle data into the appropriate cells.
Because this parallelization requires two reads of the particle data instead of one, I expected that the new version running with one thread might be slower than the old single-threaded code, especially at high particles-per-cell where the malloc overhead is lower. However, that does not appear to be the case in benchmarks. The reads are probably cheap compared to the writes, and it's possible that there are vectorization opportunities, e.g. in the computation of the cell index, that are present in this version that were not present in the old version.
Old:
New:
New, parallel:
These are direct benchmarks of
gridlink_double()(i.e. calling directly with C code; mean of 10 runs) on a 12 core, 2 socket Intel machine. Nosort_on_z, since that part is unchanged between the old and new versions. The thread scaling is far from perfect, but it's a nice little speedup. I expect this will be most helpful in latency-sensitive applications like HOD evaluations. The benchmark data is uniform random, but I've spot-checked with clustered data, which looks similar (and may be an even bigger win, because no more reallocs).The tests all pass, and I've checked with and without OpenMP and
copy_particles. Could you take a look and let me know what you think?