c++ - N体算法 : why is this slower in parallel?

Question

我编写了一个示例程序来模仿我正在处理的数据结构类型。也就是说，我有 n 对象，我需要在每个可能的对之间迭代一次并执行(对称)计算。此操作涉及将数据写入两对。在串行中，这将采用这样的循环形式

for(int i = 0; i < N-1; ++i)
   for(int j = i + 1; j < N; ++j)
      ...

但是，在 Internet 上并没有花费太多搜索就找到了“缓存不经意的并行实现”，我将其写下来并复制如下。我在这里链接了一篇详细描述该算法的帖子(使用 Intel TBB)。

https://software.intel.com/en-us/blogs/2010/07/01/n-bodies-a-parallel-tbb-solution-parallel-code-balanced-recursive-parallelism-with-parallel_invoke

我尝试使用 OpenMP 任务来执行相同的操作，但它总是比串行任务运行得慢(只是在没有 -fopenmp 的情况下进行编译)。我用 g++ -Wall -std=c++11 -O3 test.cpp -o test 编译它。使用或不使用 -O3 时都观察到同样的情况；连续剧总是更快。

要添加更多信息，在我的实际应用程序中，通常有几百到几千个元素(下面示例中的变量 n)我需要在此循环成对时尚多次。数百万次。我在下面的尝试试图模拟它(尽管我只尝试循环 10-100k 次)。

我使用 time ./test 非常粗略地计时，只是因为差异太大了。是的，我知道我的例子写得不好，而且我在我的例子中包括了创建 vector 所需的时间。但是串行时间给了我大约 30 秒，并行时间超过一分钟，所以我认为我还不需要做任何更严格的事情。

我的问题是:为什么连续剧做得更好？我在 OpenMP 中做错了什么吗？我该如何正确纠正我的错误？我是否滥用了任务？我感觉递归任务与此有关，我尝试将“OMP_THREAD_LIMIT”设置为 4，但并没有产生实质性差异。有没有更好的方法使用 OpenMP 来实现这一点？

注意:我的问题是专门询问如何修复这个特定 实现，以便它可以并行正常工作。尽管如果有人知道此问题的替代解决方案及其在 OpenMP 中的正确实现，我也对此持开放态度。

提前致谢。

#include <vector>
#include <iostream>

std::vector<std::vector<double>> testme;

void rect(int i0, int i1, int j0, int j1)
{
    int di = i1 - j0;
    int dj = j1 - j0;
    constexpr int threshold = 16;
    if(di > threshold && dj > threshold)
    {
        int im = i0 + di/2;
        int jm = j0 + dj/2;
        #pragma omp task
        { rect(i0, im, j0, jm); }
        rect(im, i1, jm, j1);
        #pragma omp taskwait

        #pragma omp task 
        { rect(i0, im, jm, j1); }
        rect(im, i1, j0, jm);
        #pragma omp taskwait
    }
    else
    {
        for(int i = i0; i < i1; ++i)
            for(int j = j0; j < j1; ++j) {
                testme[i][j] = 1;
                testme[j][i] = 1;
            }

    }
}

void triangle(int n0, int n1)
{
        int dn = n1 - n0;
        if(dn > 1)
        {
            int nm = n0 + dn/2;
            #pragma omp task
            { triangle(n0, nm); }
            triangle(nm, n1);
            #pragma omp taskwait

            rect(n0, nm, nm, n1);
        }
}


void calc_force(int nbodies)
{
    #pragma omp parallel num_threads(4)
    #pragma omp single
    triangle(0, nbodies);
}

int main()
{
    int n = 400;
    for(int i = 0; i < n; ++i)
    {
        std::vector<double> temp;
        for(int j = 0; j < n; ++j)
            temp.push_back(0);

        testme.push_back(temp);
    }

    // Repeat a bunch of times.
    for(int i = 0; i < 100000; ++i)
        calc_force(n);

    return 0;
}   

Answer 1

您当前的 OMP 任务实现似乎完全正确地应用了三角形划分方案。似乎由于分解的递归性质，当前代码只是创建了太多的子任务，调用递归三角程序直到达到 dn = 1 的条件(在树的底部)。粒度太高了。这会使您的程序负担创建和完成任务的通信要求，而任务创建的好处却越来越少；因此超过了并行性的好处。我会尝试在某个大于 1 的 dn 值(我猜大约是 15)时切断递归三角任务调用，并让最后一个(最低)任务串行执行。

线程限制只会限制事件线程的数量，但不会限制递归调用或进行的任务的数量。我会尝试一个任务 if 或者在你的三角形实现中添加一个 else。

像这样:

#include <vector>
#include <iostream>

std::vector<std::vector<double>> testme;

void rect(int i0, int i1, int j0, int j1)
{
int di = i1 - j0;
int dj = j1 - j0;
constexpr int threshold = 64;
if(di > threshold && dj > threshold)
{
    int im = i0 + di/2;
    int jm = j0 + dj/2;
    #pragma omp task
    { rect(i0, im, j0, jm); }
    rect(im, i1, jm, j1);
    #pragma omp taskwait

    #pragma omp task 
    { rect(i0, im, jm, j1); }
    rect(im, i1, j0, jm);
    #pragma omp taskwait
}
else
{
 // #pragma omp parallel for collapse(2)  (was not implimented during testing)
    for(int i = i0; i < i1; ++i)
        for(int j = j0; j < j1; ++j) {
            testme[i][j] = 1;
            testme[j][i] = 1;
        }
    }
}

void triangle(int n0, int n1)
{
    int dn = n1 - n0;
    if(dn > 1)
    {
        int nm = n0 + dn/2;
        #pragma omp task if(nm > 50 )

        { triangle(n0, nm); }
        triangle(nm, n1);
       #pragma omp taskwait

       rect(n0, nm, nm, n1);
    }
}


void calc_force(int nbodies)
{
#pragma omp parallel num_threads(4)
#pragma omp single
triangle(0, nbodies);
}

int main()
{
int n = 400;
for(int i = 0; i < n; ++i)
{
    std::vector<double> temp;
    for(int j = 0; j < n; ++j)
        temp.push_back(0);

    testme.push_back(temp);
}

// Repeat a bunch of times.
for(int i = 0; i < 100000; ++i)
    calc_force(n);

return 0;
}  

NOTE: It also might very well be the case that this implementation only shows speed up at scale, where the task overhead is outweighed by the calculation intensity of your program.