使用OpenMP并行化Ewald求和傅里叶部分的外层循环问题

我正在用OpenMP实现多处理器并行，最近尝试对Ewald求和的傅里叶部分（PME_Fourier_optimized函数）做并行化。对比并行和串行版本后发现，e_kewald能量和f_kewald力的结果在小数点后两位出现误差，而且每次运行结果都不一样，怀疑是数据竞争导致的。把外层k循环前的OpenMP并行指令去掉后，结果就和串行版本完全一致了。

我试过用#pragma omp critical解决，但不确定这是不是最优方案；后来根据建议改用fx_kewald这类1D数组结合reduction来提升并行效率，但还是不确定外层循环的OpenMP使用是否有错误，希望得到技术帮助。

参数说明 #

• ParticleN：系统带电粒子数（整数）
• kcount：整数
• qi：prop[i][0]（双精度浮点数，粒子i的电荷）
• Ext_L[0/1/2]：系统x/y/z维度（双精度浮点数）
• **r_pos：ParticleN×3数组，存储粒子坐标
• **PME_mvec_ksq：kcount×4数组
• **f_kewald：ParticleN×3数组，存储计算得到的力

实现代码 #

double PME_Fourier_optimized(int kcount, double **r_pos, double **PME_mvec_ksq, double **prop, double **f_kewald, double *Ext_L){    
    // Auto alpha determination //
    double TRTF = 50*5.5;
    double box_VOL = 8*Ext_L[0]*Ext_L[1]*Ext_L[2]; // calculate the volume of the box //
    alpha = pow((ParticleN*M_PI*M_PI*M_PI*TRTF)/(box_VOL*box_VOL),1./6.); // Ewald cutoff parameter //

    double GAMMA = -0.25/(alpha*alpha);
    double recip = 2*M_PI*ONE_PI_EPS0/(8*Ext_L[0]*Ext_L[1]*Ext_L[2]);
    double e_kewald=0.0;
    
    double fx_kewald[ParticleN];
    double fy_kewald[ParticleN];
    double fz_kewald[ParticleN];

    #pragma omp parallel for
    for (int i=0;i<ParticleN;i++){
        fx_kewald[i] = 0.0;
        fy_kewald[i] = 0.0;
        fz_kewald[i] = 0.0;
        for (int j=0;j<dimen;j++){
            f_kewald[i][j] = 0.0; 
        }
    }
    
    #pragma omp parallel for reduction(+:e_kewald,fx_kewald[:ParticleN],fy_kewald[:ParticleN],fz_kewald[:ParticleN])
    for(int k=0; k<kcount; k++){
        
        double ak_cos=0.0;
        double ak_sin=0.0;
        
        double mx = PME_mvec_ksq[k][0]; 
        double my = PME_mvec_ksq[k][1];
        double mz = PME_mvec_ksq[k][2];
        double ksq = PME_mvec_ksq[k][3];

        #pragma omp parallel for reduction(+:ak_sin,ak_cos)
        for(int i=0;i<ParticleN;i++){
            double qi = prop[i][0]; // charge of particle i //
            double kdotr = mx*r_pos[i][0] + my*r_pos[i][1] + mz*r_pos[i][2];

            ak_sin -= qi*sin(kdotr);
            ak_cos += qi*cos(kdotr);
        }
        
        double a = ak_cos;
        double b = ak_sin;

        double akak = (a*a + b*b);
        double tmp = recip * exp(GAMMA*ksq)/ksq;
        
        #pragma omp parallel for reduction (+:fx_kewald[:ParticleN],fy_kewald[:ParticleN],fz_kewald[:ParticleN])
        for(int i=0;i<ParticleN;i++){
            double qi = prop[i][0];
            double kdotr = mx*r_pos[i][0] + my*r_pos[i][1] + mz*r_pos[i][2];
            double tmp2 = 2*tmp*qi*(sin(kdotr) * a + cos(kdotr) * b);
            
            //Edit3: Following @Laci advice to use 1D array to store the forces and using the reduction to gain a huge acceleration compared to only using omp critical (See the comment section for more details).

            fx_kewald[i] += tmp2 * mx; 
            fy_kewald[i] += tmp2 * my; 
            fz_kewald[i] += tmp2 * mz;
        }
        
        e_kewald += tmp * akak; // add the energies for each k //
    }
    
    //Edit3: Finally storing the calculated values in the 2D array //
    #pragma omp parallel for
    for(int i=0;i<ParticleN;i++){
        f_kewald[i][0] = fx_kewald[i];
        f_kewald[i][1] = fy_kewald[i];
        f_kewald[i][2] = fz_kewald[i];
    }

    //printf("PME KEwald: %lf\n", e_kewald);
    return e_kewald; 
}

2D数组分配函数 #

double **alloc_2d_double(int rows, int cols){
    int i=0;
    double *data = (double *)malloc(rows*cols*sizeof(double));
    double **array= (double **)malloc(rows*sizeof(double *));
    for (i=0; i<rows; i++)
        array[i] = &(data[cols*i]);

    return array;
} 

问题根源分析 #

1. 嵌套并行冲突：外层k循环用#pragma omp parallel for创建线程池后，内层i循环又用#pragma omp parallel for重新创建线程，导致嵌套并行。这会引发线程调度混乱，内层reduction操作和外层reduction的线程上下文冲突，最终导致结果不一致。
2. 全局变量风险：代码中的alpha是全局变量，并行区域内多个线程可能同时读写该变量，引发数据竞争。
3. 数组reduction兼容性：fx_kewald[:ParticleN]这类数组reduction需要OpenMP 4.5及以上版本支持，若编译器版本过低会导致未定义行为。

修复方案 #

1. 移除嵌套并行，复用外层线程

将内层循环的#pragma omp parallel for改为#pragma omp for，让外层并行区域的线程直接处理内层任务，避免重复创建线程：

#pragma omp parallel for reduction(+:e_kewald,fx_kewald[:ParticleN],fy_kewald[:ParticleN],fz_kewald[:ParticleN])
for(int k=0; k<kcount; k++){
    // ... 省略中间代码 ...

    // 内层循环改用omp for，复用外层线程
    #pragma omp for reduction(+:ak_sin,ak_cos)
    for(int i=0;i<ParticleN;i++){
        // ... 计算ak_sin/ak_cos ...
    }

    // ... 省略中间代码 ...

    #pragma omp for reduction (+:fx_kewald[:ParticleN],fy_kewald[:ParticleN],fz_kewald[:ParticleN])
    for(int i=0;i<ParticleN;i++){
        // ... 计算力分量 ...
    }
}

2. 把全局变量alpha改为局部变量

将alpha声明为函数内的局部变量，避免线程间的读写冲突：

// 原来的全局alpha改为局部变量
double alpha = pow((ParticleN*M_PI*M_PI*M_PI*TRTF)/(box_VOL*box_VOL),1./6.);

3. 确保编译器支持OpenMP 4.5

如果使用GCC，需添加编译参数-fopenmp -std=c99或更高标准；若编译器版本过低，建议升级到GCC 5+、Clang 3.7+版本。

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