以对数时间平行减少

给定n部分和，可以将log2并行步骤中的所有部分和相加。 例如，假设有八个线程具有八个部分和： s0, s1, s2, s3, s4, s5, s6, s7 。 这可以在log2(8) = 3连续步骤中减少;

 thread0 thread1 thread2 thread4 s0 += s1 s2 += s3 s4 += s5 s6 +=s7 s0 += s2 s4 += s6 s0 += s4 

我想用OpenMP做这个，但我不想使用OpenMP的reduction条款。 我想出了一个解决方案，但我认为可以使用OpenMP的task子句找到更好的解决方案。

这比标量加法更通用。 让我选择一个更有用的案例：数组缩减（请参见此处 ， 此处以及有关数组缩减的更多信息）。

假设我想在数组a上进行数组缩减。 下面是一些代码，它们为每个线程并行填充私有数组。

 int bins = 20; int a[bins]; int **at; // array of pointers to arrays for(int i = 0; i<bins; i++) a[i] = 0; #pragma omp parallel { #pragma omp single at = (int**)malloc(sizeof *at * omp_get_num_threads()); at[omp_get_thread_num()] = (int*)malloc(sizeof **at * bins); int a_private[bins]; //arbitrary function to fill the arrays for each thread for(int i = 0; i<bins; i++) at[omp_get_thread_num()][i] = i + omp_get_thread_num(); } 

此时我有一个指向每个线程数组的指针数组。 现在我想将所有这些数组加在一起，并将最终总和写入a 。 这是我提出的解决方案。

 #pragma omp parallel { int n = omp_get_num_threads(); for(int m=1; n>1; m*=2) { int c = n%2; n/=2; #pragma omp for for(int i = 0; i<n; i++) { int *p1 = at[2*i*m], *p2 = at[2*i*m+m]; for(int j = 0; j<bins; j++) p1[j] += p2[j]; } n+=c; } #pragma omp single memcpy(a, at[0], sizeof *a*bins); free(at[omp_get_thread_num()]); #pragma omp single free(at); } 

让我试着解释一下这段代码的作用。 我们假设有八个线程。 让我们定义+=运算符来表示对数组求和。 例如s0 += s1是

 for(int i=0; i<bins; i++) s0[i] += s1[i] 

那么这段代码就行了

 n thread0 thread1 thread2 thread4 4 s0 += s1 s2 += s3 s4 += s5 s6 +=s7 2 s0 += s2 s4 += s6 1 s0 += s4 

但是这个代码并不像我想的那样理想。

一个问题是存在一些需要所有线程同步的隐式障碍。 这些障碍不是必要的。 第一个障碍是在填充arrays和进行缩减之间。 第二个障碍在#pragma omp for减少声明。 但我不能使用此方法的nowait子句来移除障碍。

另一个问题是有几个线程不需要使用。 例如，有八个线程。 减少的第一步只需要四个线程，第二步需要两个线程，最后一步只需要一个线程。 但是，这种方法将涉及减少中的所有八个线程。 虽然，其他线程无论如何都没有做多少，应该直接进入障碍并等待，这可能不是什么大问题。

我的直觉是使用omp task子句可以找到更好的方法。 不幸的是，我对task条款没什么经验，到目前为止我所做的所有努力都比我现在失败的要好。

有人可以提出一个更好的解决方案来使用例如OpenMP的task子句来减少对数时间吗？

---

我找到了一种解决障碍问题的方法。 这会异步减少。 唯一剩下的问题是它仍然将不参与减少的线程放入繁忙的循环中。 这个方法使用类似堆栈的东西将关键部分中的指针推送到堆栈（但从不弹出它们）（这是关键部分没有隐式障碍的关键之一。堆栈是串行操作但是并行减少。

这是一个有效的例子。

 #include  #include  #include  #include  void foo6() { int nthreads = 13; omp_set_num_threads(nthreads); int bins= 21; int a[bins]; int **at; int m = 0; int nsums = 0; for(int i = 0; i<bins; i++) a[i] = 0; #pragma omp parallel { int n = omp_get_num_threads(); int ithread = omp_get_thread_num(); #pragma omp single at = (int**)malloc(sizeof *at * n * 2); int* a_private = (int*)malloc(sizeof *a_private * bins); //arbitrary fill function for(int i = 0; i<bins; i++) a_private[i] = i + omp_get_thread_num(); #pragma omp critical (stack_section) at[nsums++] = a_private; while(nsums1) p1 = at[m], p2 = at[m+1], m +=2, pop = 1; if(pop) { for(int i = 0; i<bins; i++) p1[i] += p2[i]; #pragma omp critical (stack_section) at[nsums++] = p1; } } #pragma omp barrier #pragma omp single memcpy(a, at[2*n-2], sizeof **at *bins); free(a_private); #pragma omp single free(at); } for(int i = 0; i<bins; i++) printf("%d ", a[i]); puts(""); for(int i = 0; i<bins; i++) printf("%d ", (nthreads-1)*nthreads/2 +nthreads*i); puts(""); } int main(void) { foo6(); } 

我觉得使用不会在忙循环中不使用线程的任务可能会找到更好的方法。

Answer 1

实际上，使用递归的分而治之的方法实现干净利落的任务非常简单。 这几乎是教科书代码。

 void operation(int* p1, int* p2, size_t bins) { for (int i = 0; i < bins; i++) p1[i] += p2[i]; } void reduce(int** arrs, size_t bins, int begin, int end) { assert(begin < end); if (end - begin == 1) { return; } int pivot = (begin + end) / 2; /* Moving the termination condition here will avoid very short tasks, * but make the code less nice. */ #pragma omp task reduce(arrs, bins, begin, pivot); #pragma omp task reduce(arrs, bins, pivot, end); #pragma omp taskwait /* now begin and pivot contain the partial sums. */ operation(arrs[begin], arrs[pivot], bins); } /* call this within a parallel region */ #pragma omp single reduce(at, bins, 0, n); 

据我所知，没有不必要的同步，并且关键部分没有奇怪的轮询。 它也可以自然地使用不同于您的排名数据的数据大小。 我发现它非常干净，易于理解。 所以我确实认为这比你的两个解决方案都要好 。

但是让我们来看看它在实践中的表现*。 为此，我们可以使用Score-p和Vampir ：

* _{bins=10000所以减少实际上需要一点时间。 在没有涡轮增压的24核Haswell系统上执行。 gcc 4.8.4， -O3 。 我在实际执行周围添加了一些缓冲区来隐藏初始化/后处理}

该图显示了应用程序中任何线程在水平时间轴上发生的情况。 从上到下的树实现：

1. omp for loop
2. omp critical的任务。
3. omp task

这很好地展示了具体实现如何实际执行。 现在看来for循环实际上是最快的，尽管有不必要的同步。 但是，这种性能分析仍存在许多缺陷。 例如，我没有固定线程。 在实践中，NUMA（非统一内存访问）非常重要：核心是否在它自己的套接字的自己的缓存/内存中有这些数据？ 这是任务解决方案变得不确定的地方。 在简单比较中不考虑重复之间非常显着的差异。

如果还原操作在运行时变为可变，那么任务解决方案将变得比同步的for循环更好。

critical解决方案有一些有趣的方面，被动线程不会持续等待，因此它们更有可能消耗CPU资源。 这可能对性能有害，例如在turbo模式的情况下。

请记住，通过避免立即返回的产生任务， task解决方案具有更多优化潜力。 这些解决方案的执行方式也高度依赖于特定的OpenMP运行时。 英特尔的运行时似乎对任务更糟糕。

我的建议是：

• 以最佳算法复杂度实施最易维护的解决方案
• 测量代码的哪些部分对运行时实际重要
• 根据实际测量结果分析瓶颈是什么。 根据我的经验，它更多的是关于NUMA和调度，而不是一些不必要的障碍。
• 根据您的实际测量值执行微观优化

线性解决方案

以下是此问题中线性proccess_data_v1的时间轴。

OpenMP 4减少

所以我想到了OpenMP减少。 棘手的部分似乎是从循环内的数组中获取数据而没有副本。 我用NULL初始化worker数组，只是第一次移动指针：

 void meta_op(int** pp1, int* p2, size_t bins) { if (*pp1 == NULL) { *pp1 = p2; return; } operation(*pp1, p2, bins); } // ... // declare before parallel region as global int* awork = NULL; #pragma omp declare reduction(merge : int* : meta_op(&omp_out, omp_in, 100000)) initializer (omp_priv=NULL) #pragma omp for reduction(merge : awork) for (int t = 0; t < n; t++) { meta_op(&awork, at[t], bins); } 

令人惊讶的是，这看起来并不太好：

_{top是icc 16.0.2 ，bottom是gcc 5.3.0 ，都是-O3 。}

两者似乎都实现了序列化的缩减。 我试着研究一下gcc / libgomp ，但对我来说并不是很明显发生了什么。 从中间代码/反汇编，它们似乎将最终合并包装在GOMP_atomic_start / end - 这似乎是一个全局互斥锁。 类似地， icc在kmpc_critical中包含对operation的kmpc_critical 。 我认为进行昂贵的自定义减少操作没有太多优化。 传统的减少可以通过硬件支持的primefaces操作来完成。

注意每个operation如何更快，因为输入是在本地缓存的，但由于序列化，它总体上更慢。 由于差异很大，这不是一个完美的比较，早期的截图是使用不同的gcc版本。 但趋势很明显，我也有关于缓存效果的数据。