Summing over matrix columns in parallel using prange

[Andrew Hearin]

Hi,

I have written some simple code to loop over the rows of a matrix, compute the sum of the columns of each row, and return the result. The functionality should be identical to np.sum(matrix, axis=1), except I am trying to use openMP to parallelize the outer loop over the rows. My code is included below. 

The error I get is "Cannot read reduction variable in loop body". I recognize this as an openMP-specific error, and I have run a sanity check on this by removing the call to prange and the code operates as expected. I'm sure the error must be very basic since the code is so simple, but I do not see the cause when trying the parallel version. Many thanks in advance for helping with what is probably a very silly mistake. 

-----

from cython.parallel import prange, parallel

import numpy as np

cimport numpy as cnp

cimport cython

@cython.boundscheck(False)

@cython.wraparound(False)

@cython.nonecheck(False)

def parallel_sum_over_columns(cnp.float64_t[:, :] matrix):

    cdef int m = matrix.shape[0]

    cdef int n = matrix.shape[1]

    cdef int i, j

    cdef cnp.float64_t accumulator = 0.

    cdef cnp.float64_t[:] result = np.zeros(m, dtype='float64')

    with nogil, parallel(num_threads=8):

        for i in prange(m, schedule='dynamic'):

            accumulator = 0.

            for j in range(n):

                accumulator += matrix[i, j]

            result[i] = accumulator

    return np.array(result)

[Hanno Klemm]

Hi,

I think, your accumulator variable needs to be an array of size m. Otherwise, OpenMP doesn't know how to deal with a scalar variable that is defined outside of the parallel part, and then changed to different values in the different threads. 

Hanno

--

---
You received this message because you are subscribed to the Google Groups "cython-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cython-users...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Stefan Behnel]

Andrew Hearin schrieb am 02.01.2017 um 22:10:
> The error I get is "*Cannot read reduction variable in loop body*".

>
> for i in prange(m, schedule='dynamic'):
>
> accumulator = 0.
> for j in range(n):
> accumulator += matrix[i, j]
>
> result[i] = accumulator

You need to spell this out as

accumulator = accumulator + matrix[i, j]

Otherwise, Cython will assume that "accumulator" is a shared reduction
variable inside of the prange() loop.

http://docs.cython.org/en/latest/src/userguide/parallelism.html#cython.parallel.prange

Feels like we should add an example to the docs here.

Stefan

[Andrew Hearin]

Many thanks for the clear resolution to the problem. Modifying the line accumulator = accumulator + matrix[i, j]  does indeed solve the issue. 

This seems actually a bit tangential to cython and really purely openMP-related, so I'm not sure whether cython docs modifications are warranted in the short term, though it would certainly be great for the community in general to have a few minimal non-trivial worked examples to use as algorithm patterns for parallel cython development. Either way, thanks again!

Cheers,

Andrew