c++ - cython でのループの並列化

Question

ここの何人かのメンバーの助けを借りて、Python で実行され、2 つの巨大な値np.arraysを入力として受け取る関数を評価するコードを確立しました。

並列で実行されるベクトル化されたバージョンは、依然として非常に時間がかかり、シリアル fortran で記述された参照プログラムよりも約 50 倍遅くなります...

代わりに、OpenMP または MPI を使用して並列化できる cython ループを使用したいと考えています。C++ での考え方は次のようになります。

#pragma omp parallel for
for (i=0;i<np1;i++){
  for (i=0;i<np2;i++){
    double dist = sph(coord1_particle1,coord1_particle2,coord2_particle1,coord2_particle2)
    int bin=binning_function(dist)
    hist_array[bin]++
  }
}

どんなアイデアでも大歓迎です。Python のバージョンは次のとおりです。

#a is an array containing two coordinates of two objects
def dist_vec(a): # a like [[array1,array2,array2,array2],[],[]...]
  return sph(a[0],a[1],a[2],a[3]) # sph operates on coordinates

def vec_chunk(array_ab, bins) :
    dist = dist_vec(array_ab)
    hist, _ = np.histogram(dist, bins=bins)
    return hist


def mp_dist(array_a,array_b, d, bins): #d chunks AND processes
    def worker(array_ab, out_q):
        """ push result in queue """
        outdict = vec_chunk(array_ab, bins)
        out_q.put(outdict)
    # Each process will get 'chunksize' nums and a queue to put his out
    out_q = mp.Queue()
    a = np.swapaxes(array_a, 0 ,1)
    b = np.swapaxes(array_b, 0 ,1)
    array_size_a=len(array_a)-(len(array_a)%d)
    array_size_b=len(array_b)-(len(array_b)%d)
    a_chunk = array_size_a / d
    b_chunk = array_size_b / d
    procs = []
    '''prepare arrays for mp'''
    array_ab = np.empty((4, a_chunk, b_chunk))
    for j in xrange(d):
      for k in xrange(d):
        array_ab[[0, 1]] = a[:, a_chunk * j:a_chunk * (j + 1), None]
        array_ab[[2, 3]] = b[:, None, b_chunk * k:b_chunk * (k + 1)]
        p= mp.Process(target=worker, args=(array_ab, out_q))
        p.start()
        procs.append(p)
    for pro in procs:
      pro.join()
    # Collect all results into a single result dict. 
    resultarray = np.empty(len(bins)-1)
    for i in range(d):
        resultarray+=out_q.get() 
        #resultdict.update(out_q.get())
    return resultarray

bins = np.logspace(-3,1, num=25) #prepare x-axis for histogram
start_time = time()
hist_data = mp_dist(DATA,sim,10,bins)
print 'Total Time Elaspsed: ', time() - start_time