3

私は、2Dおよび3Dで地球内の変形のシミュレーションを行うために使用する、パーティクルインセルコードの並列化に取り組んでいます。コードのいくつかのルーチンは、OpenMPを使用して並列化するのが簡単で、非常にうまくスケーリングできます。ただし、粒子からグリッドセルへの補間を処理するコードの重要な部分で問題が発生しました。パーティクルは、反復ごとに動き回っています（速度フィールドに従って）。多くの計算は、通常の変形しないグリッドで実行するのに最も効率的です。したがって、各反復には、「ランダムに」分散された粒子からグリッドセルへの通信が含まれます。

この問題は、次の簡略化された1Dコードで説明できます。

//EXPLANATION OF VARIABLES (all previously allocated and initialized, 1D arrays)
//double *markerval; // Size Nm. Particle values. Are to be interpolated to the grid
//double *grid; // Size Ng=Nm/100 Grid values. 
//uint *markerpos; // Size Nm. Position of particles relative to grid (each particle
// knows what grid cell it belongs to) possible values are 0,1,...Ng-1

//#pragma omp parallel for schedule(static) private(e)
for (e=0; e<Nm; e++) {
    //#pragma omp atomic
    grid[markerpos[e]]+=markerval[e];
}

最悪のシナリオでは、パーティクルの位置はランダムですが、通常、メモリ内で互いに隣接しているパーティクルは、空間内でも、したがってグリッドメモリ内でも互いに隣接しています。

この手順を効率的に並列化するにはどうすればよいですか？複数の粒子が同じグリッドセルにマッピングされるため、上記のループが直接並列化されている場合、競合状態とキャッシュスワッピングの可能性はゼロではありません。更新をアトミックにすることで競合状態を防ぎますが、コードはシーケンシャルの場合よりもはるかに遅くなります。

また、各スレッドのグリッド値のプライベートコピーを作成し、後でそれらを合計しようとしました。ただし、これにはおそらく使用するコードに多すぎるメモリが必要であり、この例では、スレッドの数に応じて十分に拡張できませんでした（理由は不明です）。

3番目のオプションは、グリッドからパーティクルにマップしてから、パーティクルインデックスの代わりにグリッドインデックスをループすることです。ただし、これは非常に複雑で、コードを大幅に変更する必要があります。また、計算コストもかかるソートルーチンを使用する必要があるため、どれだけ役立つかはわかりません。

誰かがこれまたは同様の問題の経験を持っていますか？

4

2 に答える 2

2

オプションは、反復をスレッドに手動でマップすることです。

#pragma omp parallel shared(Nm,Ng,markerval,markerpos,grid)
{
  int nthreads = omp_get_num_threads();
  int rank     = omp_get_thread_num();
  int factor   = Ng/nthreads;

  for (int e = 0; e < Nm; e++) {
    int pos = markerpos[e];
    if ( (pos/factor)%nthreads == rank )
      grid[pos]+=markerval[e];
  }
}

いくつかのコメント:

  1. ループの反復はforスレッド間で共有されません。代わりに、各スレッドがすべての反復を行います。
  2. forループ内の条件は、どのスレッドが配列の位置posを更新するかを決定します。gridこの場所は 1 つのスレッドのみに属しているため、atomic保護は必要ありません。
  3. 式(pos/factor)%nthreadsは単なるヒューリスティックです。pos範囲内の値を返す の関数は0,...,nthreads-1、実際には、最終結果の有効性を損なうことなく、この式に置き換えることができます (より良いショットがあれば、自由に変更してください)。この機能を適切に選択しないと、負荷分散の問題が発生する可能性があることに注意してください。
于 2012-11-07T21:24:07.723 に答える
1

分子動力学アルゴリズムも OpenMP で並列化しました。まず、アルゴリズムのボトルネックを分析する必要があります(たとえば、メモリ バウンドや CPU バウンドなど)。このようにして、どこを改善すべきかがわかります。

当初、私の MD はメモリにバインドされていたので2x、データ レイアウトを構造体の配列 (AOS) から配列の構造体 (SOA) に変更するだけで速度が向上します (空間的局所性による)。また、RAM にのみ収まる入力に対して、ブロッキング手法を適用しました。元のアルゴリズムは、次のようにすべての粒子間の力のペアを計算しました。

for(int particleI = 0; i < SIZE ; i++)
 for(int particleJ = 0; j < SIZE; j++)
     calculate_force_between(i,j);

基本的にブロック法では、力の計算を粒子のブロックに集約します。たとえば、最初の 10 個の粒子間ですべての力の対を計算し、次に次の 10 個までというように計算します。

このブロック手法を使用すると、時間的局所性のより良い使用が促進されます。これは、このアプローチを使用すると、同じ粒子に対してより短い時間でより多くの計算を実現できるためです。したがって、アクセスしようとしている値がキャッシュに存在しなくなる可能性が減少します。

MD CPU がバインドされたので、 を使用して改善を試みることができmulti-threadsますが、まず、次のことを行う必要があります。

  1. アルゴリズムが実行時間のほとんどを費やしている場所を確認します。
  2. 並列で実行できるタスクを見つけて、その 粒度を決定します(並列化が正当化されるかどうかを確認するため)。
  3. Load Balance、スレッド間の作業の適切な負荷分散を保証します。
  4. 同期の使用を最小限に抑えます。

負荷分散の問題が原因で、MD のスケーリングに問題がありました。一部のスレッドは、他のスレッドより多くの作業を行っていました。ソリューション？

openMP から動的な forを試すことができます。OpenMP では、スレッドに割り当てる作業のチャンクを指定できることに注意してください。ただし、チャンクの定義には注意が必要です。dynamic forでは、チャンクが小さすぎると同期のオーバーヘッドが発生し、大きすぎると負荷分散の問題が発生する可能性があります。

また、同期のオーバーヘッドにも問題がありました。私はクリティカルを使用していましたが、アルゴリズムはスケーリングしませんでした。私はその重要性を、粒子ごとに 1 つずつ、より細かい同期、つまりロックに置き換えました。そのアプローチでいくつかの改善がありました。

(同期のオーバーヘッドに対処するための) 最後のアプローチとして、データの冗長性を使用します。各粒子はその作業を行い、結果をプライベートな一時データ構造に保存しました。最終的に、すべてのスレッドが値を減らしました。すべてのバージョンの中で、これが最高の結果をもたらしました。

CPU でかなりの高速化を達成できましたが、GPU バージョンで達成したものとは比べ物になりません。

あなたが提供した情報を使用して、私は次のようなことをします：

omp_lock_t locks [grid_size]; // create an array of locks
int g;
#pragma omp parallel for schedule(static)
for (e=0; e<Nm; e++)
{
    g = markerpos[e];

    omp_set_lock(&locks[g]);
    grid[g]+=markerval[e];
    omp_unset_lock(&locks[g]);
}

から、問題を理解したのはatomic、複数のスレッドが同じグリップ位置に同時にアクセスしないようにするために使用する必要があるということです。考えられる解決策として、ロックの配列を作成し、スレッドがグリッドの 1 つの位置にアクセスする必要があるたびに、その位置に関連付けられたロックを要求して取得することができます。別の解決策は次のとおりです。

double grid_thread[grid_size][N_threads]; // each thread have a grid
// initialize the grid_threads to zeros

#pragma omp parallel
{
    int idT = omp_get_thread_num();
    int sum;
    #pragma omp parallel for schedule(static)
    for (e=0; e<Nm; e++)
        grid_thread[markerpos[e]][idT]+=markerval[e]; // each thread compute in their 
                                                     // position
    for(int j = 0; j <Nm; j++)
    { 
        sum = 0;
        #pragma omp for reduction(+:sum) 
        for (i = 0; i < idT; i++)                   // Store the result from all
           sum += grid_thread[j][i];                // threads for grid position j

         #pragma barrier                            // Ensure mutual exclusion

         #pragma master
         grid[j] +=sum;                             // thread master save the result  
                                                    // original grid
         #pragma barrier                            // Ensure mutual exclusion
      }
   }
}
于 2012-11-07T23:42:19.993 に答える