OpenMP和Cython缺乏加速和错误的结果

Question

我正在尝试使用OpenMP加速用Cython编写的一段简单代码.这是一个双循环,对于输入数组中的每个位置,在每个参考点添加一个数量.这是代码的主要部分:

cimport cython  
import numpy as np
cimport numpy as np
cimport openmp
DTYPE = np.double
ctypedef np.double_t DTYPE_t

cdef extern from "math.h" nogil :
  DTYPE_t sqrt(DTYPE_t)


@cython.cdivision(True)
@cython.boundscheck(False)
def summation(np.ndarray[DTYPE_t,ndim=2] pos, np.ndarray[DTYPE_t,ndim=1] weights, 
              np.ndarray[DTYPE_t, ndim=2] points, int num_threads = 0):

    from cython.parallel cimport prange, parallel, threadid

    if num_threads <= 0 : 
        num_threads = openmp.omp_get_num_procs()

    if num_threads > openmp.omp_get_num_procs() : 
        num_threads = openmp.omp_get_num_procs()

    openmp.omp_set_num_threads(num_threads)

    cdef unsigned int nips = len(points)
    cdef np.ndarray[DTYPE_t, ndim=1] sum_array = np.zeros(nips, dtype = np.float64)
    cdef np.ndarray[DTYPE_t, ndim=2] sum_array3d = np.zeros((nips,3), dtype = np.float64)

    cdef unsigned int n = len(weights)

    cdef unsigned int pi, i, id
    cdef double dx, dy, dz, dr, weight_i, xi,yi,zi

    print 'num_threads = ', openmp.omp_get_num_threads()

    for i in prange(n,nogil=True,schedule='static'):
        weight_i = weights[i]
        xi = pos[i,0]
        yi = pos[i,1]
        zi = pos[i,2]
        for pi in range(nips) :
            dx = points[pi,0] - xi
            dy = points[pi,1] - yi
            dz = points[pi,2] - zi
            dr = 1.0/sqrt(dx*dx + dy*dy + dz*dz)
            sum_array[pi] += weight_i * dr
            sum_array3d[pi,0] += weight_i * dx
            sum_array3d[pi,1] += weight_i * dy
            sum_array3d[pi,2] += weight_i * dz


    return sum_array, sum_array3d

我已将其与相关的测试和设置文件放在一个要点(https://gist.github.com/rokroskar/6ed1bfc1a5f8f9c183a6)

出现两个问题:

首先,在当前配置中,我无法获得任何加速.代码在多个内核上运行,但时间显示没有优势.

其次,结果根据表示存在竞争条件的核心数而不同.原地资金是否应该最终减少？或者是有趣的事情发生,因为它是一个嵌套的循环？我认为prange每个线程中的所有内容都是单独执行的.

如果我颠倒循环的顺序,这两个都会消失 - 但是因为我的外循环现在是结构化的方式,所有数据读取完成的地方,如果我反转它们,数组遍历num_thread次,这是浪费的.我也尝试将整个嵌套循环放在一个with parallel():块中并显式使用线程局部缓冲区,但无法使其工作.

很显然,我遗漏了一些关于OpenMP应该如何工作的基本知识(虽然这可能是Cython特有的？)所以我会很感激提示!

Answer 1

您是否尝试过切换两个循环，以免多个线程从同一位置读取和写入？我非常确定这些循环不会自动提升为 OpenMP 缩减，并且诸如“sum_array3d[pi,0] += Weight_i * dx”之类的增量不是原子的。

此外，由于计算相对简单，Cython 可能会太过分，您可能会使用Parakeet或Numba来代替。

默认情况下，Parakeet 将使用 OpenMP 并行执行推导式。您必须重写代码，使其看起来像：

@parakeet.jit
def summation(pos, weights, points):
  n_points = len(points)
  n_weights = len(weights)
  sum_array3d = np.zeros((n_points,3))
  def compute(i):
    pxi = points[i, 0]
    pyi = points[i, 1]
    pzi = points[i, 2]
    total = 0.0
    for j in xrange(n_weights):
      weight_j = weights[j]
      xj = pos[j,0]
      yj = pos[j,1]
      zj = pos[j,2]
      dx = pxi - pos[j, 0]
      dy = pyi - pos[j, 1]
      dz = pzi - pos[j, 2]
      dr = 1.0/np.sqrt(dx*dx + dy*dy + dz*dz)
      total += weight_j * dr
      sum_array3d[i,0] += weight_j * dx
      sum_array3d[i,1] += weight_j * dy
      sum_array3d[i,2] += weight_j * dz
    return total 
  sum_array = np.array([compute(i) for i in xrange(n_points)])
  return sum_array, sum_array3d

至于 Numba，我不确定prange结构是否已经进入免费版本。

编辑：额头拍打，抱歉，我错过了您考虑切换循环的问题部分。