Parallel prefix sum for OpenMP backend

[Slaedr]

Simple algorithm with O(n log n) work and O(log n) depth or span.

Edit: It was decided to change the algorithm to one with O(n) work and O(p) depth (p = number of processors). Even though it's less parallel in a sense, it now has less memory accesses.

[upsj]

Great idea, we can really use this. I would suggest a slightly modified algorithm though, which would reduce the amount of synchronization and memory access by adding a small sequential part:

• Each thread gets a block of data (like right now)
• Each thread computes the local prefix sum + total block sum
• Each thread writes the total block sum into total_sums[thread_id]
• One thread computes the prefix-sum over total_sums
• Each thread adds total_sums[thread_id] to each of its local entries.

This would mean 2n reads and 2n writes in parallel plus p reads and writes sequentially for p threads

[Slaedr]

Great idea, we can really use this. I would suggest a slightly modified algorithm though, which would reduce the amount of synchronization and memory access by adding a small sequential part:

* Each thread gets a block of data (like right now)

* Each thread computes the local prefix sum + total block sum

* Each thread writes the total block sum into `total_sums[thread_id]`

* One thread computes the prefix-sum over total_sums

* Each thread adds `total_sums[thread_id]` to each of its local entries.

This would mean 2n reads and 2n writes in parallel plus p reads and writes sequentially for p threads

Thanks for taking a look at it. I had actually thought of this first and it's simpler to implement. The advantage is the total work would be O(n) instead of O(n log n), but the disadvantage is the max depth would be O(p). In a weak scaling scenario, which to me is the most important for exascale etc., the max depth would be O(n). The extra memory allocation also becomes O(n). The current approach might be slower for some small problem sizes, but might be better suited for Ginkgo, no? The current approach also avoids dynamic memory allocation altogether.

[upsj]

Can you clarify what you mean by depth? It sounds like you are implying a size limit for this approach? If necessary, the prefix sum on total_sums can also be parallelized if the number of cores is too large (even though this sounds unlikely on-node). Finally, we only need O(p) additional memory, so this should really be negligible.

[upsj]

Suggested change

	if (num_entries <= 1) {
	return;
	}
	if (num_entries <= 1) {
	counts[0] = 0;
	return;
	}

[Slaedr]

Can you clarify what you mean by depth? It sounds like you are implying a size limit for this approach? If necessary, the prefix sum on total_sums can also be parallelized if the number of cores is too large (even though this sounds unlikely on-node). Finally, we only need O(p) additional memory, so this should really be negligible.

By depth I mean the number of sequential steps in the kernel. You make a good point - because this is CPUs we're talking about, the number of cores should typically be very small. In the largest cases maybe it's up to 100 or so but seldom more than that.

But by having an extra array in main memory, would we not be introducing more latency? Right now, what you would be storing in the global total_sums array, I just store in a local register. I think pretty much all of my accesses are cache-friendly. (There would be a few cache-line conflicts which would happen in both implementations, but that's possibly avoidable in the current implementation by moving the boundaries of work chunks a bit.) I wonder if the extra cost of 1 logarithm and log n number of power operations are worth avoiding the O(p) uncached main memory accesses.

[upsj]

The reads and writes into total_sums should be negligible, since after the small prefix sum computation, each thread only needs to read a single element from total_sums, so it can be kept in a register, too. If you look at line 109, if I see this correctly, each element is read and written O(log p) times instead of 1 time in the two-level approach.

[Slaedr]

format!

[codecov]

Codecov Report

Merging #703 (89914a7) into develop (14da01e) will increase coverage by 0.00%.
The diff coverage is 100.00%.

@@           Coverage Diff            @@
##           develop     #703   +/-   ##
========================================
  Coverage    92.44%   92.45%           
========================================
  Files          363      362    -1     
  Lines        26923    26920    -3     
========================================
- Hits         24890    24888    -2     
+ Misses        2033     2032    -1     

Impacted Files	Coverage Δ
omp/test/components/prefix_sum.cpp	100.00% <100.00%> (ø)
include/ginkgo/core/base/range_accessors.hpp	78.94% <0.00%> (+5.26%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 14da01e...89914a7. Read the comment docs.

[upsj]

LGTM! Only small suggestions

[upsj]

I don't think we should restrict ourselves here: There are some cases (conversion between empty matrices, ...) where we call prefix_sum on 0 or 1-sized arrays, the expected output would not be a failure, but writing 0 or nothing to the output array.

[fritzgoebel]

LGTM

[pratikvn]

LGTM!

[sonarcloud]

Kudos, SonarCloud Quality Gate passed!

0 Bugs
0 Vulnerabilities
0 Security Hotspots
2 Code Smells

No Coverage information
0.0% Duplication