C++ OpenMP并行for循环使其比单线程慢得多。

Question

我对C中的OpenMP是个新手。我用它来在函数中并行我的for循环，但事实证明，与单线程的情况相比，它大大减慢了我的for循环。例如，对于每个点(光环)，For循环可以在10s左右完成，但使用OpenMP需要几分钟。

在这个函数中，我试图计算每个点(光晕)的许多壳的密度，通过计数壳内的粒子，然后将它们写到一个数组中。我要计算的粒子有512^3个，大约有200个点(光晕)。我想拆分不同线程的点(光环)，以使其更快。

#define ArrayAccess2D_n2(a, n1, n2, i1, i2) (a)[ i2+n2*i1 ]


void halo_shell_rho(float boxsize, float *halo_pos, float *halo_R, int halo_number,\ 
int halo_start, int halo_end, float *par_pos, long long par_number,\ 
int shell_bins, float rmax_fac, float *out_shell_den){

    float temp;

    long long iter_sfs=0, iter_sfc=0, iter_ufs=0, iter_ufc=0;
    int dim=3;

    float par_posx, par_posy, par_posz, dist;
    float halo_posx, halo_posy, halo_posz, halo_rad;
    int i=0, ini_j=0, vol_j=0;
    int a=0, b=0;
    long long k=0;

    #pragma omp parallel for private(i, ini_j, vol_j, a, b, k)
    for(i=halo_start; i<=halo_end; i++){
            printf("halo %d\n", i);
            float count[shell_bins];
            float volume[shell_bins];

            for(ini_j=0; ini_j<shell_bins; ini_j++){
                    count[ini_j] = 0;
                    volume[ini_j] = 0; }

            halo_posx = ArrayAccess2D_n2(halo_pos, dim, halo_number, 0, i);
            halo_posy = ArrayAccess2D_n2(halo_pos, dim, halo_number, 1, i);
            halo_posz = ArrayAccess2D_n2(halo_pos, dim, halo_number, 2, i);
            halo_rad = halo_R[i];

            for(vol_j=0; vol_j<shell_bins; vol_j++){

                    volume[vol_j] = shell_volume((vol_j+1)*halo_rad*rmax_fac/(shell_bins*1000), vol_j*halo_rad*rmax_fac/(shell_bins*1000)); }

            for(k=0; k<par_number; k++){

                    par_posx = ArrayAccess2D_n2(par_pos, par_number, dim, k, 0);
                    par_posy = ArrayAccess2D_n2(par_pos, par_number, dim, k, 1);
                    par_posz = ArrayAccess2D_n2(par_pos, par_number, dim, k, 2);

                    dist = pb_distance(boxsize*1000, halo_posx, halo_posy, halo_posz, par_posx, par_posy, par_posz); //1000 for boxsize in Mpc

                    if(dist <= 2*rmax_fac*halo_rad){

                            for(a=0; a<shell_bins; a++){

                                    if((dist <= halo_rad*(a+1)*rmax_fac/shell_bins) && (dist >= halo_rad*a*rmax_fac/shell_bins)){

                                            count[a] += 1; }
                            }
                    }
            }

            for(b=0; b<shell_bins; b++){

            out_shell_den[(i-halo_start+0*(1+halo_end-halo_start))*shell_bins+b] = count[b]/volume[b]; 
            //out_shell_den has shape (2, halo_number, shell_bins), 0 for edge, 1 for density
            out_shell_den[(i-halo_start+1*(1+halo_end-halo_start))*shell_bins+b] = (2*b+1)*rmax_fac/(shell_bins*2);
            }
    }

}

有人能帮我解决这个问题吗？我知道这是一个经常被问到的问题，但我没有从其他帖子中找到任何解决方案。如果有帮助，我在一个有32个线程的集群上运行它。

谢谢!

Answer 1

感谢@DavidSchwartz和@tim18。

像halo_rad和par_posx这样的变量是在并行之前声明的，这意味着它们被隐式假设为公共的。所以它变慢了，因为所有的线程都在争夺它们的使用权。解决这个问题的一种方法是将所有变量添加到private()中。但我认为更好的方法是在并行中声明变量，如下所示：

void halo_shell_rho(float boxsize, float *halo_pos, float *halo_R, int halo_number, int halo_start, int halo_end, float *par_pos, long long par_number, int shell_bins, float rmax_fac, float *out_shell_den){

    int dim=3;
    int i=0, ini_j=0, vol_j=0, a=0, b=0;
    long long k=0;

    #pragma omp parallel for private(i, ini_j, vol_j, a, b, k)
    for(i=halo_start; i<=halo_end; i++){
            printf("halo %d\n", i);

            float halo_posx, halo_posy, halo_posz, halo_rad;
            float count[shell_bins];
            float volume[shell_bins];

            for(ini_j=0; ini_j<shell_bins; ini_j++){
                    count[ini_j] = 0;
                    volume[ini_j] = 0; }

            halo_posx = ArrayAccess2D_n2(halo_pos, dim, halo_number, 0, i);
            halo_posy = ArrayAccess2D_n2(halo_pos, dim, halo_number, 1, i);
            halo_posz = ArrayAccess2D_n2(halo_pos, dim, halo_number, 2, i);
            halo_rad = halo_R[i];

            for(vol_j=0; vol_j<shell_bins; vol_j++){

                    volume[vol_j] = shell_volume((vol_j+1)*halo_rad*rmax_fac/(shell_bins*1000), vol_j*halo_rad*rmax_fac/(shell_bins*1000)); }


            for(k=0; k<par_number; k++){
                    float par_posx, par_posy, par_posz, dist;

                    par_posx = ArrayAccess2D_n2(par_pos, par_number, dim, k, 0);
                    par_posy = ArrayAccess2D_n2(par_pos, par_number, dim, k, 1);
                    par_posz = ArrayAccess2D_n2(par_pos, par_number, dim, k, 2);

                    dist = pb_distance(boxsize*1000, halo_posx, halo_posy, halo_posz, par_posx, par_posy, par_posz); //1000 for boxsize in Mpc

                    if(dist <= 2*rmax_fac*halo_rad){

                            for(a=0; a<shell_bins; a++){

                                    if((dist <= halo_rad*(a+1)*rmax_fac/shell_bins) && (dist >= halo_rad*a*rmax_fac/shell_bins)){

                                            count[a] += 1; }
                            }
                    }
            }

            for(b=0; b<shell_bins; b++){

            out_shell_den[(i-halo_start+0*(1+halo_end-halo_start))*shell_bins+b] = count[b]/volume[b]; //out_shell_den has shape (2, halo_number, shell_bins), 0 for edge, 1 for density
            out_shell_den[(i-halo_start+1*(1+halo_end-halo_start))*shell_bins+b] = (2*b+1)*rmax_fac/(shell_bins*2);
            }
    }
}