C OpenMP parallel forループはシングルスレッドよりもはるかに遅くなります

Question

私はOpenMPをC言語で使い慣れていますが、関数内でforループを並列化するのに使用しましたが、場合。たとえば、forループは各ポイント（ハロー）に対して約10秒で実行できますが、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のような変数は、パラレルの前に宣言されています。つまり、暗黙的にパブリックとみなされます。したがって、すべてのスレッドがスレッドを使用する権利を超えて戦っているため、速度が低下します。これを解決する1つの方法は、すべての変数を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); 
      } 
    } 
}