我正在学习如何使用OpenMP使代码使用多个处理器。最近，我尝试使用OpenMP使埃瓦尔德求和傅立叶部分并行。下面是名为PME_fourier_oprimized的函数，我尝试使用OpenMP杂注使其并行，其中，

ParticleN = Number of cahrged particles in the System (an integer number)
kcount = an integer number
qi = prop[i][0] (a double float) which gives charge of a particle i
Ext_L[0] = x dimension of the system (a double float number)
Ext_L[1] = y dimension of the system (a double float number)
Ext_L[1] = z dimension of the system (a double float number)
**r_pos = ParticleNx3 array (r_pos[i][0], r_pos[i][1], and r_pos[i][2] = x, y , and z postions of a particle i)
**PME_mev_ksq = kcountx4 array
**f_kewald = ParticleNx3 array which stores the forces

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; 
}

我将这个代码输出与我的串行代码版本进行了比较，注意到e_kewald能量和f_kewald力在小数点后2位有错误。有趣的是，每次运行它都得到不同的结果（可能是数据竞争的情况）。
当我删除外层k循环之前的第一个pragma命令（#pragma omp parallel for private（k，mx，my，mz，ksq，qi，kdotr，tmp2）reduction（+：e_kewald））时，它开始完美地工作并与串行代码结果完美匹配。我是否做错了什么？我也不确定我使用pargma的方式是否正确。任何帮助都是非常感谢的。
Edit1：在f_kewald数组解决问题之前使用#pragma omp critical之后。但是，不确定这是否是最好的方法。
我的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;
}

我仍然想知道f_kewald [i][0]，f_kewald [i][1]，f_kewald [i][2]是否可以在约简下使用。
Edit2：声明循环中的变量，并在omp临界区中放置大括号。
微粒数量~5000 - 10000千个~2000
Edit3：使用带约简的1D数组（fx_kewald，fy_kewald，fz_kewald）来提高并行化速度，而不是使用带2D数组的opm临界区。（感谢@Laci下面的评论）