How to do parallel array `.sum`?

[fzyzcjy]

Hi thanks for the lib! It seems that array.sum() is a bit slow. So I wonder whether there is a way to parallelly sum it? I have lookat at parallel but it seems only to support map, and does not support things like .fold().

[adamreichold]

Generally, when folding in parallel, Rayon folds within a given task and then reduces the accumulators to a single value, c.f. https://docs.rs/rayon/1.5.1/rayon/iter/trait.ParallelIterator.html#method.fold

To apply this to e.g. summing the elements of an array, you need to access the underlying parallel iterations, for example

let sum = array.view()
    .into_par_iter() 
    .fold(
        || 0.0,
        |mut acc, val| {
            acc += val;
            acc
        },
    )
    .reduce(
        || 0.0,
        |mut lhs, rhs| {
            lhs += rhs;
            lhs
        }
    );

I am not sure, but since someone asked me this recently, I think it could make sense for this library to expose this common case as a dedicated function to improve discoverability of the idiom? Maybe something like

fn par_fold<ID, F, R, T>(self, identity: ID, fold: F, reduce: R) -> T
where
    ID: Fn() -> T + Send + Sync + Clone,
    F: Fn(T, Self::Item) -> T + Send + Sync,
    R: Fn(T, T) -> T + Send + Sync,
    T: Send
{
    self.into_par_iter().fold(identity.clone(), fold).reduce(identity, reduce)
}

?

(One thing that is important to keep in mind in numerical applications is that the positions of splitting and reduction are non-deterministic so that for exact results, the operations need to be associative which floating point operation are not. So for floating point operations the above will imply larger errors.)

[bluss]

The guideline I have is that sum (which element type?!) for floats is memory bound, lots of memory accesses compared with the number of arith ops. There should not be anything to gain from thread parallelism. I would look first towards data layout and simd. Also, implement something and get to know how it really behaves.

[bluss]

@adamreichold If the reduction for example would split the array in half every time and sum each half, then you expect to reduce the rounding error compared to a linear sum.

We don't really like the float += loop since it doesn't vectorize when written plainly unfortunately. ndarray's .sum() tries to do the vectorization part better, if it can.

[adamreichold]

The guideline I have is that sum (which element type?!) for floats is memory bound, lots of memory accesses compared with the number of arith ops. There should not be anything to gain from thread parallelism.

I think this makes sense for an element-wise sum but the fold-reduce operation could have a very different granularity operating on e.g. generalized rows or blocks which are processed sequentially within the closure.

If the reduction for example would split the array in half every time and sum each half, then you expect to reduce the rounding error compared to a linear sum.

Agreed, geometric splitting is better for summing. I think the main point I was trying and failing to make is that the grouping of accumulators is not deterministic and hence an analysis of precision must also consider the worst-case splitting for any given operation.

We don't really like the float += loop since it doesn't vectorize when written plainly unfortunately. ndarray's .sum() tries to do the vectorization part better, if it can.

If the library provided a par_sum method based on par_fold, it could use the internal numeric_util module. But I don't think par_sum is particularly useful, it is rather the fold-reduce idiom which I think would be good to make more discoverable as it can be applied at a more useful coarser granularity than element-wise.

[fzyzcjy]

Thanks for the information! I will try them.