c - 对数时间并行减少

Question

给定 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 找到更好的解决方案。条款。

这比标量加法更通用。让我选择一个更有用的案例:数组缩减(有关数组缩减的更多信息，请参阅 here 、 here 和 here )。

假设我想对数组进行数组缩减 a .这是一些为每个线程并行填充私有(private)数组的代码。

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

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

一个问题是有一些隐式障碍需要所有线程同步。这些障碍不应该是必要的。第一个障碍是填充数组和减少之间。第二个障碍在 #pragma omp for减持声明。但我不能使用 nowait子句用这个方法来消除障碍。

另一个问题是有几个线程不需要使用。例如有八个线程。还原的第一步只需要四个线程，第二步两个线程，最后一步只需要一个线程。但是，此方法将涉及所有八个线程的减少。尽管如此，其他线程无论如何都不会做太多事情，应该直接进入屏障并等待，所以这可能不是什么大问题。

我的直觉是使用 omp task 可以找到更好的方法。条款。不幸的是，我对 task 几乎没有经验。条款和我迄今为止的所有努力都比我现在失败的减少更好。

有人可以建议一个更好的解决方案来减少对数时间，例如使用OpenMP的 task条款？

我找到了一种解决障碍问题的方法。这会异步减少。唯一剩下的问题是它仍然将不参与减少的线程放入繁忙的循环中。此方法使用类似堆栈的东西将指针插入临界区中的堆栈(但从不弹出它们)(这是 critical sections don't have implicit barriers 的键之一。堆栈串行操作但并行减少。

这是一个工作示例。

#include <stdio.h>
#include <omp.h>
#include <stdlib.h>
#include <string.h>

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(nsums<2*n-2) {
            int *p1, *p2;
            char pop = 0;
            #pragma omp critical (stack_section)
            if((nsums-m)>1) 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

实际上，使用递归分而治之的方法通过任务干净利落地实现这一点非常简单。这差不多textbook代码。

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 .我在实际执行周围添加了一些缓冲区以隐藏初始化/后处理



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

omp for循环

omp critical一种任务。

omp task

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

如果减少操作在运行时变得可变，那么任务解决方案将变得比同步 for 循环更好。
critical解决方案有一些有趣的方面，被动线程不会持续等待，因此它们更有可能消耗 CPU 资源。这可能对性能不利，例如在涡轮模式的情况下。

请记住 task通过避免立即返回的生成任务，解决方案具有更大的优化潜力。这些解决方案的性能还很大程度上取决于特定的 OpenMP 运行时。英特尔的运行时似乎对任务的处理要差得多。

我的建议是:

用最优算法实现最易维护的解决方案
复杂性

衡量代码的哪些部分对运行时真正重要

根据实际测量分析什么是瓶颈。根据我的经验，它更多地是关于 NUMA 和调度，而不是一些不必要的障碍。

根据您的实际测量进行微优化

线性解决方案

这是线性 proccess_data_v1 的时间表来自 this question .



OpenMP 4 缩减

所以我想到了 OpenMP 减少。棘手的部分似乎是从 at 获取数据。循环内的数组，没有副本。我确实用 NULL 初始化了工作数组并在第一次简单地移动指针:

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);
        }

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



顶部是 icc 16.0.2 , 底部是 gcc 5.3.0 , 都与 -O3 .

两者似乎都实现了序列化的减少。我试图查看 gcc/ libgomp ，但对我来说发生的事情并不是很明显。从中间代码/反汇编来看，他们似乎将最终合并包装在 GOMP_atomic_start 中。/ end - 这似乎是一个全局互斥锁。同样 icc将调用包装到 operation在 kmpc_critical .我想在昂贵的自定义减少操作中没有太多优化。传统的减少可以通过硬件支持的原子操作来完成。

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