Cython: Странная ошибка выполнения параллельно, но без параллели

Question

---

В Cython я пытаюсь сделать FOR LOOP параллельно. Итак, я компилирую его, используя ниже, как упоминалось везде:

ext_type = Extension("wrap_aip",
                 sources=["wrap_aip.pyx", "my_aip.c"],
                 extra_compile_args=['-fopenmp'],
                 xtra_link_args=['-fopenmp'],
                 include_dirs=[numpy.get_include()])

Код может выполняться параллельно, но результаты вычислений областей выдают несколько ошибочных значений. (Например, 86 из 1000 полигонов неправильны)
После многих проверок я только комментирую -fopenmp, как показано ниже:

ext_type = Extension("wrap_aip",
                 sources=["wrap_aip.pyx", "my_aip.c"],
                 # extra_compile_args=['-fopenmp'],
                 # extra_link_args=['-fopenmp'],
                 include_dirs=[numpy.get_include()])

С вычислительными результатами областей все в порядке, но он не может работать параллельно!!! Я использую Memoryviews и Pointer для передачи данных в блок nogil, как упомянуто в Cython-doc. Разве это не безопасно параллельно? Как передать данные в блок nogil?
И у меня есть другая версия, передающая данные по np.PyArray_DATA с помощью gil. Это всегда правильно, но не может работать параллельно. Код части FOR LOOP показан ниже:

# poly: np.array in size of N*28, means N polygons.
# wrap_aip_d.intersection: function computing the area of intersection given a pair of polygons.

def polygon_NMS(polys, scores, thresh=0.7):
scores = np.array(scores)
order = np.argsort(scores)[::-1]
cdef:
    int i = 0
    int j = 0
    int num_polys = polys.shape[0]
    float inS = -1.
    np.ndarray[float, mode="c", ndim=2] inter_areas = np.zeros((num_polys, num_polys), dtype=np.float32)
    np.ndarray[float, mode="c", ndim=1] areas = np.zeros(num_polys, dtype=np.float32)
    float [:,:] f_poly=polys.astype(np.float32)
    float * poly_i
    float * poly_j

## Important Part, nogil block ##
for i in prange(num_polys,nogil=True):
    for j in prange(i, num_polys):
        poly_i=&f_poly[i][0]
        poly_j=&f_poly[j][0]
        areas[i] = wrap_aip_d.intersection(poly_i, poly_i)
        inS = wrap_aip_d.intersection(poly_i, poly_j)
        inter_areas[i, j] = inS
        inter_areas[j, i] = inS

## Not important below! ##

keep = []
while order.size > 0:
    i = order[0]
    keep.append(i)

    ovr = inter_areas[i][order[1:]] / (areas[i] + areas[order[1:]] - inter_areas[i][order[1:]])
    inds = np.where(ovr <= thresh)[0]
    order = order[inds + 1]
return keep, areas, inter_areas

Если вы хотите проверить больше кода:
Версия 2:(Правильно все время, но не может работать параллельно)

def gil_pnms(polys, scores, thresh=0.7):
scores = np.array(scores)
order = np.argsort(scores)[::-1]
cdef:
    int i = 0
    int j = 0
    int num_polys = polys.shape[0]
    float inS = -1.
    np.ndarray[float, mode="c", ndim=2] ps = polys.astype(np.float32)
    np.ndarray[float, mode="c", ndim=2] inter_areas = np.zeros((num_polys, num_polys), dtype=np.float32)
    np.ndarray[float, mode="c", ndim=1] areas = np.zeros(num_polys, dtype=np.float32)

for i in prange(num_polys,nogil=True):
    for j in prange(i, num_polys):
        with gil:
            a = <float*> np.PyArray_DATA(ps[i])
            b = <float*> np.PyArray_DATA(ps[j])
        areas[i] = wrap_aip_d.intersection(a, a)
        inS = wrap_aip_d.intersection(a, b)
        inter_areas[i, j] = inS
        inter_areas[j, i] = inS

keep = []
while order.size > 0:
    i = order[0]
    keep.append(i)
    ovr = inter_areas[i][order[1:]] / (areas[i] + areas[order[1:]] - inter_areas[i][order[1:]])
    inds = np.where(ovr <= thresh)[0]
    order = order[inds + 1]
return keep,areas,inter_areas