对这样的工作负载使用递归算法的简单想法对我来说已经非常陌生了。然后，使用OpenMP任务对其进行并行化似乎更陌生...为什么不使用更传统的方法来解决这个问题呢？
所以我决定用我想到的几种方法给予一试。但是为了让练习变得有意义，重要的是要为计算“对称计算”做一些实际的工作，否则，仅仅在一个所有元素上迭代而不考虑对称方面肯定是最好的选择。
所以我写了一个a force()函数，根据两个物体的坐标计算与它们之间引力相互作用松散相关的东西，然后我测试了4种不同的方法来迭代粒子：
1.像你提出的那种简单的三角形方法，由于它内在的负载不平衡特性，这个方法被一个schedule(auto)子句并行化，以允许运行时库做出它认为对性能最好的决定。
1.一个“聪明”的三角形域遍历，包括在j方向上将其切成两半，以允许使用2个规则的循环。基本上，它对应于如下内容：

/|
  / |       __  __
 /  |  =>  | //   |
/___|      |//____|

1.一个简单的矩形方法，只是忽略对称性。注意，这个方法，就像你的递归方法一样，保证了对force数组的非并发访问。
1.一种线性化方法，包括预先计算用于访问三角域的索引i和j的阶，并且迭代包含这些索引的向量。
由于其中利用force[i] += fij; force[j] -= fij;方法来累积力的向量将为非并行化索引中的更新生成竞态条件（例如方法#1中的j），我创建了一个本地线程前强制数组，在进入并行区域时初始化为0，然后在线程前对这个“私有”数组进行计算。并且在退出并行区域时，个体贡献被累积在具有critical构造的全局力阵列上，这是OpenMP中阵列的典型缩减模式。
下面是完整的代码：

#include <iostream>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <omp.h>

std::vector< std::vector<double> > l_f;
std::vector<double> x, y, f;
std::vector<int> vi, vj;

int numth, tid;
#pragma omp threadprivate( tid )

double force( int i, int j ) {
    double dx = x[i] - x[j];
    double dy = y[i] - y[j];
    double dist2 = dx*dx + dy*dy;
    return dist2 * std::sqrt( dist2 );
}

void loop1( int n ) {
    #pragma omp parallel
    {
        for ( int i = 0; i < n; i++ ) {
            l_f[tid][i] = 0;
        }
        #pragma omp for schedule(auto) nowait
        for ( int i = 0; i < n-1; i++ ) {
            for ( int j = i+1; j < n; j++ ) {
                double fij = force( i, j );
                l_f[tid][i] += fij;
                l_f[tid][j] -= fij;
            }
        }
        #pragma omp critical
        for ( int i = 0; i < n; i++ ) {
            f[i] += l_f[tid][i];
        }
    }
}

void loop2( int n ) {
    int m = n/2-1+n%2;
    #pragma omp parallel
    {
        for ( int i = 0; i < n; i++ ) {
            l_f[tid][i] = 0;
        }
        #pragma omp for schedule(static) nowait
        for ( int i = 0; i < n; i++ ) { 
            for ( int j = 0; j < m; j++ ) {
                int ii, jj;
                if ( j < i ) {
                    ii = n-1-i;
                    jj = n-1-j;
                }
                else {
                    ii = i;
                    jj = j+1;
                }
                double fij = force( ii, jj );
                l_f[tid][ii] += fij;
                l_f[tid][jj] -= fij;
            }
        }
        if ( n%2 == 0 ) {
            #pragma omp for schedule(static) nowait
            for ( int i = 0; i < n/2; i++ ) {
                double fij = force( i, n/2 );
                l_f[tid][i] += fij;
                l_f[tid][n/2] -= fij;
            }
        }
        #pragma omp critical
        for ( int i = 0; i < n; i++ ) {
            f[i] += l_f[tid][i];
        }
    }
}

void loop3( int n ) {
    #pragma omp parallel for schedule(static)
    for ( int i = 0; i < n; i++ ) {
        for ( int j = 0; j < n; j++ ) {
            if ( i < j ) {
                f[i] += force( i, j );
            }
            else if ( i > j ) {
                f[i] -= force( i, j );
            }
        }
    }
}

void loop4( int n ) {
    #pragma omp parallel
    {
        for ( int i = 0; i < n; i++ ) {
            l_f[tid][i] = 0;
        }
        #pragma omp for schedule(static) nowait
        for ( int k = 0; k < vi.size(); k++ ) {
            int i = vi[k];
            int j = vj[k];
            double fij = force( i, j );
            l_f[tid][i] += fij;
            l_f[tid][j] -= fij;
        }
        #pragma omp critical
        for ( int i = 0; i < n; i++ ) {
            f[i] += l_f[tid][i];
        }
    }
}

int main( int argc, char *argv[] ) {
    if ( argc != 2 ) {
        std::cout << "need the dim as argument\n";
        return 1;
    }
    int n = std::atoi( argv[1] );

    // initialise data
    f.resize( n );
    x.resize( n );
    y.resize( n );
    for ( int i = 0; i < n; ++i ) {
        x[i] = y[i] = i;
        f[i] = 0;
    }

    // initialising linearised index vectors
    for ( int i = 0; i < n-1; i++ ) {
        for ( int j = i+1; j < n; j++ ) {
            vi.push_back( i );
            vj.push_back( j );
        }
    }
    // initialising the local forces vectors
    #pragma omp parallel
    {
        tid = omp_get_thread_num();
        #pragma master
        numth = omp_get_num_threads();
    }
    l_f.resize( numth );
    for ( int i = 0; i < numth; i++ ) {
        l_f[i].resize( n );
    }

    // testing all methods one after the other, with a warm up before each
    int lmax = 10000;
    loop1( n );
    double tbeg = omp_get_wtime();
    for ( int l = 0; l < lmax; l++ ) {
        loop1( n );
    }
    double tend = omp_get_wtime();
    std::cout << "Time for triangular loop is " << tend-tbeg << "s\n";

    loop2( n );
    tbeg = omp_get_wtime();
    for ( int l = 0; l < lmax; l++ ) {
        loop2( n );
    }
    tend = omp_get_wtime();
    std::cout << "Time for mangled rectangular loop is " << tend-tbeg << "s\n";

    loop3( n );
    tbeg = omp_get_wtime();
    for ( int l = 0; l < lmax; l++ ) {
        loop3( n );
    }
    tend = omp_get_wtime();
    std::cout << "Time for naive rectangular loop is " << tend-tbeg << "s\n";

    loop4( n );
    tbeg = omp_get_wtime();
    for ( int l = 0; l < lmax; l++ ) {
        loop4( n );
    }
    tend = omp_get_wtime();
    std::cout << "Time for linearised loop is " << tend-tbeg << "s\n";

    int ret = f[n-1];
    return ret;
}

现在，评估它们的相对性能和可伸缩性变得简单了，所有方法在第一次非计时预热迭代之后，都在循环中计时。
编制：

g++ -O3 -mtune=native -march=native -fopenmp tbf.cc -o tbf

8核IvyBridge CPU上的结果：

> OMP_NUM_THREADS=1 numactl -N 0 -m 0 ./tbf 500
Time for triangular loop is 9.21198s
Time for mangled rectangular loop is 10.1316s
Time for naive rectangular loop is 15.9408s
Time for linearised loop is 10.6449s
> OMP_NUM_THREADS=2 numactl -N 0 -m 0 ./tbf 500
Time for triangular loop is 6.84671s
Time for mangled rectangular loop is 5.13731s
Time for naive rectangular loop is 8.09542s
Time for linearised loop is 5.4654s
> OMP_NUM_THREADS=4 numactl -N 0 -m 0 ./tbf 500
Time for triangular loop is 4.03016s
Time for mangled rectangular loop is 2.90809s
Time for naive rectangular loop is 4.45373s
Time for linearised loop is 2.7733s
> OMP_NUM_THREADS=8 numactl -N 0 -m 0 ./tbf 500
Time for triangular loop is 2.31051s
Time for mangled rectangular loop is 2.05854s
Time for naive rectangular loop is 3.03463s
Time for linearised loop is 1.7106s

因此，在这种情况下，方法#4似乎是最好的选择，同时具有良好的性能和非常好的可伸缩性。请注意，由于schedule(auto)指令提供了良好的负载平衡工作，直接的三角形方法并不太差。但最终，我鼓励您使用自己的工作负载进行测试...
作为参考，您的初始代码（修改后用于计算force()，方法与其他测试完全相同，包括使用的OpenMP线程数，但不需要本地强制数组和最终归约，而是使用递归方法）如下所示：

> OMP_NUM_THREADS=1 numactl -N 0 -m 0 ./recursive 500
Time for recursive method is 9.32888s
> OMP_NUM_THREADS=2 numactl -N 0 -m 0 ./recursive 500
Time for recursive method is 9.48718s
> OMP_NUM_THREADS=4 numactl -N 0 -m 0 ./recursive 500
Time for recursive method is 10.962s
> OMP_NUM_THREADS=8 numactl -N 0 -m 0 ./recursive 500
Time for recursive method is 13.2786

您当前的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;
}

注意：也可能出现这种情况，即这种实现只显示了大规模的速度提升，其中任务开销超过了程序的计算强度。

并行代码没有发挥其潜力的原因之一是由于一个称为false sharing (Wikipedia)的问题。
解决方案是对问题进行分区，以便输出2-D矩阵中的每个缓存行（内部向量）仅由一个线程更新。这对于通过构造的三角形是成立的，但是如果im和jm不是在高速缓存行边界处对准的条目的索引，则矩形中的并行性是有问题的。如果im和/或jm指示的分区不在缓存线边界，则两个线程将写入公共缓存线，但该高速缓存线内的偏移量不同-假共享的定义。
This article by Intel对虚假共享有很好的描述，以及如何避免的建议。
https://software.intel.com/en-us/articles/avoiding-and-identifying-false-sharing-among-threads
现引述有关章节，以供参考：
伪共享是SMP系统上众所周知的性能问题，SMP系统中的每个处理器都有一个本地缓存。当不同处理器上的线程修改驻留在同一缓存线上的变量时，就会发生这种情况，如图1所示。这种情况称为伪共享，因为每个线程实际上并没有共享对同一变量的访问。访问同一变量或真共享，将需要编程同步构造来确保有序的数据访问。
下面的示例代码中以红色显示的源代码行会导致错误共享：

double sum=0.0, sum_local[NUM_THREADS];
#pragma omp parallel num_threads(NUM_THREADS)
{
 int me = omp_get_thread_num();
 sum_local[me] = 0.0;

 #pragma omp for
 for (i = 0; i < N; i++)
 sum_local[me] += x[i] * y[i];

 #pragma omp atomic
 sum += sum_local[me];
}

数组sum_local可能存在错误共享。此数组的大小取决于线程数，并且足够小，可以放入单个缓存行。并行执行时，线程会修改sum_local（源行显示为红色）的不同但相邻的元素，这会使所有处理器该高速缓存行无效。
x1c 0d1x图1.当不同处理器上的线程修改驻留在同一缓存线上的变量时，会发生假共享。这会使该高速缓存线无效，并强制更新内存以保持缓存一致性。
在图1中，线程0和1需要内存中相邻且驻留在同一缓存线上的变量。该高速缓存线被加载到CPU 0和CPU 1的缓存中（灰色箭头）。即使线程修改了不同的变量（红色和蓝色箭头），缓存线也会失效，从而强制更新内存以保持缓存一致性。
Here's a lecture介绍了使用OpenMP组织N-Body算法的各种方法的优缺点，但请注意Walter在评论中的警告，本讲座更多的是关于编程而不是物理，正如Walter在N体问题中所指出的，在该问题中，一些粒子与力的计算紧密接触，以确定物体上的净力（加速度）和积分以确定速度，然后再次确定位置，必须仔细完成-全局阶跃函数是不合适的。
http://www.cs.usask.ca/~spiteri/CMPT851/notes/nBody.pdf
特别是第18页，我将其复制在此以供参考：

应用Foster的方法因此，任务到核心的Map简化为粒子到核心的Map。假设每个步骤完成的工作大致相等，为每个核心分配大致n/p个粒子的块分区应提供良好平衡的负载。此假设适用于计算fi，j时未利用对称性的情况（t）。当利用对称性时，用于较小i的循环将比用于较大i的循环更昂贵。在这种情况下，循环划分更有效。然而，在共享存储器框架中，循环分区几乎肯定会导致比块分区更高数量的高速缓存未命中。在分布式存储器框架中，循环分区所涉及的通信开销将可能大于块分区所涉及的通信开销

基于任务的并行性的艺术在于避免订阅不足和订阅过多。这意味着在某个时候必须串行执行任务，因为并行执行由于开销而变得太慢（另请参见discussion here）。何时达到这一点取决于工作量和任务调度程序。
在rect()函数中，我们已经使用了threshold来限制任务并行执行到每边有超过threshold个元素的区域，但奇怪的是，在triangle()中没有这样做，所以我的第一个攻击方法是在该例程中也尝试类似的技术。

void triangle(int n0, int n1, const int threshold)
{
    int dn = n1 - n0;
    if(dn > threshold)
    {
        int nm = n0 + dn/2;
        #pragma omp task
        { triangle(n0, nm, threshold); }
        triangle(nm, n1, threshold);
        #pragma omp taskwait
        rect(n0, nm, nm, n1, threshold);  // pass threshold on
    } else {
        for(int i = n0; i < n1; ++i)
            for(int j = i+1; j < n1; ++j) {   // excludes self-interactions
                auto fij  = mutual_force(i,j);
                force[i] += fij;
                force[j] -= fij;
            }
    }
}

注意，我将threshold作为一个运行时变量。这允许您试验它，以查看定时对它的依赖程度。通常的依赖关系是一个长谷，具有良好的缩放性，但值太大或太小会产生不好的结果。对于一个好的任务调度器，您希望生成比线程多得多的任务，但threshold也要比1大得多，比如64-1024。
当然，这两个要求之间存在挤压：你不能有效地将小问题扩展到许多线程，强扩展有它的限制（N个操作不能在多于N个线程之间共享）。
很可能是因为你的问题太小，无法通过这种方式高效地并行化，特别是只有几百个粒子。另一种并行化策略是计算每对粒子的力两次，并使用简单的基于for循环的并行化

#pragma omp parallel for
for(int i=0; i<n; ++i) {
    force[i] = 0;
    for(int j=0; j<n; ++j)
        force[i] += mutual_force(i,j);
}

编译器会发现这很容易优化，静态并行也可以。

我猜你看到了次优的缓存使用率。现在的内存比CPU慢得多--多个数量级。
当线程A占用字节1，线程B占用字节2，线程C占用字节3时，三个CPU内核将分别需要将一个完整的缓存行放入其L1缓存中，才能使用其中的一个字节。CPU必须确保缓存是一致的，并且由于缓存行开头的数据在其他任何内容之前可用，它们将以不同的速度前进。然后可能开始颠簸公共的更高级高速缓存。
另一方面，在单线程版本中，一个CPU内核必须完成所有工作，但它可以获得最佳的可用内存访问：可预测的升序访问，并且总是使用整个高速缓存行。
怎么解决这样的问题呢？那么你已经开始做了，性能工程最重要的是什么：测量。也许您可以通过确保每个线程处理一个完整的缓存线，而另一个线程处理另一个缓存线来解决这个问题。我不知道OpenMP是否支持对线程上的工作负载分布进行细粒度控制。测量一下，看看是否有帮助。