在 C 中使用 OpenMP 进行外循环矢量化的问题

问：

我正在研究如何使用 OpenMP 使代码使用多个处理器。最近，我尝试使用 OpenMP 使我的 Ewald 求和傅里叶零件并行。下面是我尝试使用 OpenMP 编译指示并行PME_fourier_oprimized的函数，其中，

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

我将此代码输出与我的序列号版本进行了比较，发现小数点后 2 位后的能量和f_kewald力e_kewald存在错误。有趣的是，每次运行它时我都会得到不同的结果（也许是数据竞争的情况）。

当我在外部 k 循环之前删除第一个 pragma 命令 （ #pragma omp parallel for private（k， mx， my， mz， ksq， qi， kdotr， tmp2） reduction（+：e_kewald） ） 时，它开始完美工作并与串行代码结果完美匹配。我做错了什么吗？我也不确定我使用 pargma 的方式是否正确。任何帮助都非常感谢。

编辑1：在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] 是否可以在还原下使用。

编辑2：在循环中声明变量，并将大括号放在omp关键部分。

颗粒N ~ 5000-10000 kcount ~ 2000

编辑3：使用1D数组（fx_kewald，fy_kewald，fz_kewald）来获得增强的并行化速度，而不是使用2D数组的opm关键部分。（感谢下面@Laci的评论）