使用OpenMP静态调度并行化双层for循环的问题求助

问题描述 #

尝试用OpenMP静态调度并行化N体问题的加速度计算函数，但结果不正确，怀疑是accelerations[i]的伪共享问题导致。

原始串行代码 #

void computeAccelerations(){
int i,j;
for(i=0;i<bodies;i++){
    accelerations[i].x = 0; accelerations[i].y = 0; accelerations[i].z = 0;
    for(j=0;j<bodies;j++){
        if(i!=j){
            //accelerations[i] = addVectors(accelerations[i],scaleVector(GravConstant*masses[j]/pow(mod(subtractVectors(positions[i],positions[j])),3),subtractVectors(positions[j],positions[i])));
            vector sij = {positions[i].x-positions[j].x,positions[i].y-positions[j].y,positions[i].z-positions[j].z};
            vector sji = {positions[j].x-positions[i].x,positions[j].y-positions[i].y,positions[j].z-positions[i].z};
            double mod = sqrt(sij.x*sij.x + sij.y*sij.y + sij.z*sij.z);
            double mod3 = mod * mod * mod;
            double s = GravConstant*masses[j]/mod3;
            vector S = {s*sji.x,s*sji.y,s*sji.z};
            accelerations[i].x+=S.x;accelerations[i].y+=S.y;accelerations[i].z+=S.z;
        }
    }
}

并行尝试代码 #

void computeAccelerations_static(int num_of_threads){
int i,j;
#pragma omp parallel for num_threads(num_of_threads) schedule(static)
for(i=0;i<bodies;i++){
    accelerations[i].x = 0; accelerations[i].y = 0; accelerations[i].z = 0;
    for(j=0;j<bodies;j++){
        if(i!=j){
            //accelerations[i] = addVectors(accelerations[i],scaleVector(GravConstant*masses[j]/pow(mod(subtractVectors(positions[i],positions[j])),3),subtractVectors(positions[j],positions[i])));
            vector sij = {positions[i].x-positions[j].x,positions[i].y-positions[j].y,positions[i].z-positions[j].z};
            vector sji = {positions[j].x-positions[i].x,positions[j].y-positions[i].y,positions[j].z-positions[i].z};
            double mod = sqrt(sij.x*sij.x + sij.y*sij.y + sij.z*sij.z);
            double mod3 = mod * mod * mod;
            double s = GravConstant*masses[j]/mod3;
            vector S = {s*sji.x,s*sji.y,s*sji.z};
            accelerations[i].x+=S.x;accelerations[i].y+=S.y;accelerations[i].z+=S.z;
        }
    }
}

问题分析与修复方案 #

伪共享的原因 #

accelerations数组中相邻的vector元素大概率处于同一CPU缓存行。当多个线程同时修改不同i对应的accelerations[i]时，缓存行会频繁因缓存一致性协议失效，既拖慢性能，也可能导致计算结果异常。

核心修复：线程私有临时变量 #

每个线程先在私有临时变量中完成当前i的加速度计算，最后一次性写入全局数组，避免频繁修改全局内存引发的伪共享：

void computeAccelerations_static(int num_of_threads){
    int i,j;
    #pragma omp parallel for num_threads(num_of_threads) schedule(static) private(j)
    for(i=0;i<bodies;i++){
        // 线程私有临时变量，存储当前i的加速度计算结果
        vector acc_temp = {0.0, 0.0, 0.0};
        for(j=0;j<bodies;j++){
            if(i!=j){
                vector sij = {positions[i].x-positions[j].x, positions[i].y-positions[j].y, positions[i].z-positions[j].z};
                double mod = sqrt(sij.x*sij.x + sij.y*sij.y + sij.z*sij.z);
                double mod3 = mod * mod * mod;
                double s = GravConstant*masses[j]/mod3;
                // 简化计算：sji = -sij，直接用-s*sij替代s*sji
                acc_temp.x += -s * sij.x;
                acc_temp.y += -s * sij.y;
                acc_temp.z += -s * sij.z;
            }
        }
        // 最后一次性写入全局数组，减少全局内存访问次数
        accelerations[i] = acc_temp;
    }
}

额外优化建议 #

• 简化计算逻辑：利用sji = -sij的关系，省去sji和S的定义，减少内存开销与计算步骤。
• 显式私有变量声明：在OpenMP指令中添加private(j)，确保每个线程拥有独立的循环变量j，避免潜在的线程冲突。
• 缓存行对齐：如果仍存在性能瓶颈，可给vector结构体添加缓存行对齐属性，确保每个accelerations元素独占一个缓存行（以64字节缓存行为例）：

  typedef struct {
      double x;
      double y;
      double z;
  } vector __attribute__((aligned(64)));
  

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