openmp並列ループの2倍の削減*

Question

配列（double **）に対する縮小（合計）のためのOpenMPの使用に問題があります。コードは次のとおりです。

#include <stdlib.h> 
#include <iostream> 
#include <omp.h> 

using namespace std ; 

double my_function(const double * x) { return 1 ; } 

int main() 
{ 
    const int nb_mom_x = 5 ; 
    const int nb_mom_y = 5 ; 
    const int mesh_nb_pts_x = 1000 ; 
    const int mesh_nb_pts_y = 1000 ; 
    double PDF_moments[nb_mom_x][nb_mom_y] ; 
    double tmp_coordinates[2] ; 

    // Initialisation : 
    for (int k_ordre_mom_x = 0 ; k_ordre_mom_x < nb_mom_x ; k_ordre_mom_x++) 
     for(int k_ordre_mom_y = 0 ; k_ordre_mom_y < nb_mom_y ; k_ordre_mom_y++) 
      PDF_moments[k_ordre_mom_x][k_ordre_mom_y] = 0.0 ; 

    #pragma omp parallel for reduction(+:PDF_moments) 
    for (int k_mesh_x = 0 ; k_mesh_x < mesh_nb_pts_x ; k_mesh_x++) 
    { 
     // Abscissa of the point in x: 
     tmp_coordinates[0] = (double) k_mesh_x/mesh_nb_pts_x ; 

     // Evaluation of the powers of x: 
     double tmp_pow_x[nb_mom_x] ; 
     tmp_pow_x[0] = 1.0 ; 
     for (int k_ordre_mom_x = 1 ; k_ordre_mom_x < nb_mom_x ; k_ordre_mom_x++) 
      tmp_pow_x[k_ordre_mom_x] = tmp_pow_x[k_ordre_mom_x-1] * tmp_coordinates[0] ; 

     for (int k_mesh_y = 0 ; k_mesh_y < mesh_nb_pts_y ; k_mesh_y++) 
     { 
      // Abscissa of the point in y: 
      tmp_coordinates[1] = (double) k_mesh_y/mesh_nb_pts_y ; 

      // Evaluation of the powers of y: 
      double tmp_pow_y[nb_mom_y] ; 
      tmp_pow_y[0] = 1.0 ; 
      for (int k_ordre_mom_y = 1 ; k_ordre_mom_y < nb_mom_x ; k_ordre_mom_y++) 
       tmp_pow_y[k_ordre_mom_y] = tmp_pow_y[k_ordre_mom_y-1] * tmp_coordinates[0] ; 

      // Evaluation of the function: 
      double tmp_function = my_function(tmp_coordinates) ; 

      for(int k_ordre_mom_y = 0 ; k_ordre_mom_y < nb_mom_y ; k_ordre_mom_y++) 
       for (int k_ordre_mom_x = 0 ; k_ordre_mom_x < nb_mom_x ; k_ordre_mom_x++) 
        PDF_moments[k_ordre_mom_x][k_ordre_mom_y] += tmp_function * tmp_pow_x[k_ordre_mom_x] * tmp_pow_y[k_ordre_mom_y] ; 

     } 
    } 

    for (int k_ordre_mom_x = 0 ; k_ordre_mom_x < nb_mom_x ; k_ordre_mom_x++) 
     for(int k_ordre_mom_y = 0 ; k_ordre_mom_y < nb_mom_y ; k_ordre_mom_y++) 
      PDF_moments[k_ordre_mom_x][k_ordre_mom_y] *= (1.0/mesh_nb_pts_x)*(1.0/mesh_nb_pts_y) ;  


    // Display of the moments of the function: 
    printf("\n+-------+") ; 
    for (int k_ordre_mom_y = 0 ; k_ordre_mom_y < nb_mom_y ; k_ordre_mom_y++) printf("------------+") ; 
    printf("\n| Ordre |") ; 
    for (int k_ordre_mom_y = 0 ; k_ordre_mom_y < nb_mom_y ; k_ordre_mom_y++) printf(" % 9.1d |",k_ordre_mom_y) ; 
    printf("\n+-------+") ; 
    for (int k_ordre_mom_y = 0 ; k_ordre_mom_y < nb_mom_y ; k_ordre_mom_y++) printf("------------+") ; 
    for (int k_ordre_mom_x = 0 ; k_ordre_mom_x < nb_mom_x ; k_ordre_mom_x++) 
    { 
     printf("\n| % 5.1d |",k_ordre_mom_x) ; 
     for(int k_ordre_mom_y = 0 ; k_ordre_mom_y < nb_mom_y ; k_ordre_mom_y++) 
     { 
      printf(" % 1.3e |" ,PDF_moments[k_ordre_mom_x][k_ordre_mom_y]) ; 
     } 
    } 
    printf("\n+-------+") ; 
    for (int k_ordre_mom_y = 0 ; k_ordre_mom_y < nb_mom_y ; k_ordre_mom_y++) printf("------------+") ; 

    printf("\n"); 
    return 0 ; 
} 

問題は、私は結果が悪事あるMac上でg++-7 -fopenmp main.cpp -O0 -o out.exeで、このコードをコンパイルするときは、次のです。しかし、私がこのコードのコンパイル後に#pragma omp parallel for reduction(+:PDF_moments)という行にコメントすると、結果は正しいです。

私の推測では、コマンドのオプションは：#pragma omp parallel for reduction(+:PDF_moments)を修正する必要がありますが、それを行う方法はわかりません。

Answer 1

tmp_coordinatesに競合状態があります。ベクタは、各ループ反復で独立していても、ループの外側で宣言されるため、デフォルトで共有されます。代わりに、ループの内側に宣言します。

#pragma omp parallel for reduction(+:PDF_moments) 
for (int k_mesh_x = 0 ; k_mesh_x < mesh_nb_pts_x ; k_mesh_x++) 
{ 
    double tmp_coordinates[2] ; 

次に、スレッドごとにプライベートです。