OpenMP并行N体程序线程数超2后性能未提升问题排查

我用C++实现了时间复杂度为O(N²)的朴素N体问题，并用OpenMP尝试并行化。测试参数为1000个天体、100个时间步长，多次运行取中位数的结果如下：

• 单线程：14.06565秒
• 2线程：6.25671秒，加速比甚至超过2倍
• 3线程：6.90014秒，相比2线程无任何提升
• 4线程：7.66918秒，线程数超过4后性能也没有改善

我的CPU是AMD Ryzen 5 3600，核心数量足够支持更多线程同时运行。之前用OpenMP从未遇到过仅2线程就达到最优性能且加速比超2倍的情况，因此怀疑代码存在问题，但无法定位具体原因。我已经尝试过以下操作，性能表现完全一致：

• 调整OpenMP调度子句
• 移除可能存在伪共享的变量/数组（比如二维向量f）
• 注释掉moveBodies函数

问题似乎集中在calculateForces函数，怀疑是负载均衡问题？

代码如下：

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

using namespace std;

struct point{
    double x=0;
    double y=0;

    point(double x_, double y_) : x(x_), y(y_) {}

    point() : x(0), y(0) {}
};

double G = 6.67e-11; //gravity

#define DT 0.01
#define SIZE 9e7


//Generate random double
double fRand(double fMin, double fMax) {
    double f = (double)rand() / RAND_MAX;
    return fMin + f * (fMax - fMin);
}


//calculate total force for every pair of bodies
void calculateForces(vector<point>& p, vector<vector<point>>& f, vector<double>& m, const int& n){
    double distance;
    double magnitude;
    point direction;

    #pragma omp parallel for schedule(static, 1) //Stripes evenly distributed
    for (int i=0; i<n-1;i++){
        int id = omp_get_thread_num();
        for (int j = i+1; j < n; j++) {
            distance =  sqrt(pow((p[i].x - p[j].x),2) + pow((p[i].y - p[j].y),2));
            magnitude = (G*m[i]*m[j]) / pow(distance,2);
            direction.x = p[j].x-p[i].x; 
            direction.y = p[j].y-p[i].y;

            f[id][i].x = f[id][i].x + magnitude*direction.x/distance;

            f[id][j].x = f[id][j].x - magnitude*direction.x/distance;

            f[id][i].y = f[id][i].y + magnitude*direction.y/distance;
            
            f[id][j].y = f[id][j].y - magnitude*direction.y/distance;
        }
    }
}

//calculate new velocity and position for each body
void moveBodies(vector<point>& p, vector<vector<point>>& f, vector<double>& m, vector<point>& v, const int& n) {
    point deltav; // dv=f/m * DT
    point deltap; // dp=(v+dv/2) * DT

    #pragma omp parallel for private(deltav, deltap) schedule(static, 1)
    for (int i = 0; i < n; i++) {
        point force;
        for (int k = 0; k < omp_get_num_threads(); k++) {
            force.x += f[k][i].x; f[k][i].x = 0;
            force.y += f[k][i].y; f[k][i].y = 0;
        }
        deltav.x = force.x/m[i] * DT;
        deltav.y = force.y/m[i] * DT;
        deltap.x = (v[i].x + deltav.x/2) * DT;
        deltap.y = (v[i].y + deltav.y/2) * DT;
        v[i].x = v[i].x + deltav.x;
        v[i].y = v[i].y + deltav.y;
        p[i].x = p[i].x + deltap.x;
        p[i].y = p[i].y + deltap.y;
        force.x = force.y = 0.0; // reset force vector
    }
}

int main(int argc, char** argv)
{

    if (argc != 4) {
        cout << "Incorrect amount of arguments!" << endl;
        return 0;
    }

    srand(time(NULL));

    int gnumBodies = atoi(argv[1]);
    int numSteps = atoi(argv[2]);
    int nThreads = atoi(argv[3]);

    cout << "number of bodies: " << gnumBodies << endl;
    cout << "number of timesteps: " << numSteps << endl;

    omp_set_num_threads(nThreads);

    vector<point> p; //Position
    vector<point> v; //velocity
    vector<vector<point>> f(nThreads, vector<point>(gnumBodies)); //force
    vector<double> m; //mass

    for (int i = 0; i < gnumBodies; i++) {
        p.push_back(point(fRand(0, SIZE), fRand(0, SIZE)));
        v.push_back(point(0,0));
        m.push_back(fRand(6e22, 2e30));
    }

    auto start = omp_get_wtime();
    for (int i=0; i<numSteps; i++){
        calculateForces(p, f, m, gnumBodies);
        moveBodies(p, f, m, v, gnumBodies);
    }

    cout << "Time taken was " << (omp_get_wtime() - start)<< " s"<< endl;
}

内容的提问来源于stack exchange，提问作者EyedBread