python - 这个问题可以在 Cython 中与 OpenMP 并行实现吗？

Question

我已将一些 Cython 代码与 OpenMP 并行化。 偶尔，代码会计算出错误的结果。

我为我的问题创建了一个几乎最小的工作示例。 “几乎”，因为错误结果的频率似乎取决于代码中最微小的更改，因此，例如我保留了函数指针。

Cython代码是

#cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
# distutils: language = c++

import numpy as np

cimport cython
from cython.parallel import prange, parallel
from libcpp.vector cimport vector
cimport numpy as np

cdef inline double estimator_matheron(const double f_diff) nogil:
    return f_diff * f_diff

ctypedef double (*_estimator_func)(const double) nogil

cdef inline void normalization_matheron(
    vector[double]& variogram,
    vector[long]& counts,
    const int variogram_len
):
    cdef int i
    for i in range(variogram_len):
        if counts[i] == 0:
            counts[i] = 1
        variogram[i] /= (2. * counts[i])

ctypedef void (*_normalization_func)(vector[double]&, vector[long]&, const int)


def test(const double[:] f):
    cdef _estimator_func estimator_func = estimator_matheron
    cdef _normalization_func normalization_func = normalization_matheron

    cdef int i_max = f.shape[0] - 1
    cdef int j_max = i_max + 1

    cdef vector[double] variogram_local, variogram
    cdef vector[long] counts_local, counts

    cdef int i, j

    with nogil, parallel():
        variogram_local.resize(j_max, 0.0)
        counts_local.resize(j_max, 0)

        for i in range(i_max):
            for j in range(1, j_max-i):
                counts_local[j] += 1
                variogram_local[j] += estimator_func(f[i] - f[i+j])

    normalization_func(variogram_local, counts_local, j_max)

    return np.asarray(variogram_local)

为了测试代码，我使用了这个脚本:

import numpy as np
from cython_parallel import test

z = np.array(
    (41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3),
    dtype=np.double,
)

print(test(z))

结果应该是

[0.         0.49166667 0.7625     1.09071429 0.90166667 1.336
 0.9525     0.435      0.005      0.405     ]

这是错误结果的典型表现

[0.         0.44319444 0.75483871 1.09053571 0.90166667 1.336
 0.9525     0.435      0.005      0.405     ]

此代码主要将数字求和到向量 variogram_local 中。大多数时候，这段代码是有效的，但如果没有进行足够的统计，可能每 30 次就会产生错误的结果。如果我将 with nogil, parallel(): 行更改为 with nogil:，它总是有效。如果我根本不使用函数指针，它也总是有效，如下所示:

    with nogil, parallel():
        variogram_local.resize(j_max, 0.0)
        counts_local.resize(j_max, 0)

        for i in range(i_max):
            for j in range(1, j_max-i):
                counts_local[j] += 1
                variogram_local[j] += (f[i] - f[i+j]) * (f[i] - f[i+j])

    for j in range(j_max):
        if counts_local[j] == 0:
            counts_local[j] = 1
        variogram_local[j] /= (2. * counts_local[j])

    return np.asarray(variogram_local)

完整的代码在不同的平台上进行了测试，这些问题主要发生在 macOS 上，有 clang，例如:

https://ci.appveyor.com/project/conda-forge/staged-recipes/builds/29018878

编辑

感谢您的意见，我修改了代码并使用 num_threads=2 它可以正常工作。但是一旦 num_threads>2 我又得到了错误的结果。您是否认为，如果 Cython 对 OpenMP 的支持是完美的，我的新代码应该可以工作，还是我仍然出错？如果这是站在 Cython 一边，我想我确实会用纯 C++ 实现代码。

def test(const double[:] f):
    cdef int i_max = f.shape[0] - 1
    cdef int j_max = i_max + 1

    cdef vector[double] variogram_local, variogram
    cdef vector[long] counts_local, counts

    cdef int i, j, k

    variogram.resize(j_max, 0.0)
    counts.resize(j_max, 0)

    with nogil, parallel(num_threads=2):
        variogram_local = vector[double](j_max, 0.0)
        counts_local = vector[long)(j_max, 0)

        for i in prange(i_max):
            for j in range(1, j_max-i):
                counts_local[j] += 1
                variogram_local[j] += (f[i] - f[i+j]) * (f[i] - f[i+j])

        for k in range(j_max):
            counts[k] += counts_local[k]
            variogram[k] += variogram_local[k]

    for i in range(j_max):
        if counts[i] == 0:
            counts[i] = 1
        variogram[i] /= (2. * counts[i])

    return np.asarray(variogram)

Answer 1

与其名称相反，variogram_local 和 counts_local 实际上不是本地的。它们是共享的，所有线程并行处理它们，因此是未定义的结果。

请注意，您实际上并没有分享任何工作。只是所有线程都在做同样的事情 - 整个串行任务。

比较合理的并行版本看起来更像这样:

variogram.resize(j_max, 0.0)
counts.resize(j_max, 0)

with nogil, parallel():
    for i in range(i_max):
        for j in prange(1, j_max-i):
            counts[j] += 1
            variogram[j] += estimator_func(f[i] - f[i+j])

共享数组在外部初始化，然后线程共享内部 j 循环。因为没有两个线程会在同一个 j 上工作，所以这样做是安全的。

现在将内循环并行化可能并不理想。如果您实际上要并行化外部循环，您实际上必须创建实际的局部变量并在之后合并/减少它们。