Cython shared memory в cython.parallel.prange - блок

У меня есть функция foo, которая берет указатель на память как аргумент и записывает и читает в эту память:

cdef void foo (double *data):
   data[some_index_int] = some_value_double
   do_something_dependent_on (data)

Я выделяю data так:

cdef int N = some_int
cdef double *data = <double*> malloc (N * sizeof (double))

cdef int i
for i in cython.parallel.prange (N, nogil=True):
    foo (data)

readout (data)

Теперь мой вопрос: как разные темы рассматривают это? Я предполагаю, что память, на которую указывает data, будет разделяться всеми потоками и "одновременно" читать или записывать, находясь внутри функции foo. Это могло бы испортить все результаты, поскольку нельзя полагаться на ранее установленную дату-дату (в пределах foo)? Является ли мое предположение правильным или есть какой-то волшебный защитный пояс, реализованный в компиляторе cython?

Спасибо вам заблаговременно.

Ответ 1

Хороший способ состоит в том, чтобы основной массив был изолирован от нитей. Затем вы даете каждому потоку указатель на часть основного массива, которая должна быть вычислена потоком.

Следующий пример представляет собой реализацию матричного умножения (аналогично dot для двухмерных массивов), где:

c = a*b

parallelism здесь реализуется по строкам a. Проверьте, как указатели передаются в функцию multiply, чтобы разные потоки могли использовать одни и те же массивы.

import numpy as np
cimport numpy as np
import cython
from cython.parallel import prange

ctypedef np.double_t cDOUBLE
DOUBLE = np.float64


def mydot(np.ndarray[cDOUBLE, ndim=2] a, np.ndarray[cDOUBLE, ndim=2] b):
    cdef np.ndarray[cDOUBLE, ndim=2] c
    cdef int i, M, N, K

    c = np.zeros((a.shape[0], b.shape[1]), dtype=DOUBLE)
    M = a.shape[0]
    N = a.shape[1]
    K = b.shape[1]

    for i in prange(M, nogil=True):
        multiply(&a[i,0], &b[0,0], &c[i,0], N, K)

    return c


@cython.wraparound(False)
@cython.boundscheck(False)
@cython.nonecheck(False)
cdef void multiply(double *a, double *b, double *c, int N, int K) nogil:
    cdef int j, k
    for j in range(N):
        for k in range(K):
            c[k] += a[j]*b[k+j*K]

Чтобы проверить, вы можете использовать этот script:

import time

import numpy as np

import _stack

a = np.random.random((10000,500))
b = np.random.random((500,2000))

t = time.clock()
c = np.dot(a, b)
print('finished dot: {} s'.format(time.clock()-t))

t = time.clock()
c2 = _stack.mydot(a, b)
print('finished mydot: {} s'.format(time.clock()-t))

print 'Passed test:', np.allclose(c, c2)

Где на моем компьютере он дает:

finished dot: 0.601547366526 s
finished mydot: 2.834147917 s
Passed test: True

Если число строк a было меньше, чем число cols или число cols в b, то mydot было бы хуже, требуя лучшей проверки того, какое измерение делает parallelism.

Ответ 2

Я предполагаю, что без блокировки синхронизации чтения или записи в data потоки будут читать/записывать в ячейку памяти и перезаписывать друг друга. Вы не получите согласованных результатов без какой-либо синхронизации.

Хотя документы (http://docs.cython.org/src/userguide/parallelism.html), похоже, предполагают, что OpenMP (бэкэнд по умолчанию) автоматически создает локаторы потоков.