1

OpenMP に問題があります。並列ブロックで何かをインクリメントする場合、その式の前にアトミックを設定する必要があることを知っています。しかし、私のコードには理解できない部分があります。

ここでアトミックを使用する必要があるのはなぜですか?

#pragma omp parallel
{
  double distance, magnitude, factor, r;
  vector_t direction;
  int i, j;
#pragma omp for
  for (i = 0; i < n_body - 1; i++)
    {
      for (j = i + 1; j < n_body; j++)
        {
          r = SQR (bodies[i].position.x - bodies[j].position.x) + SQR (bodies[i].position.y - bodies[j].position.y);
          // avoid numerical instabilities
          if (r < EPSILON)
            {
              // this is not how nature works :-)
              r += EPSILON;
            }
          distance = sqrt (r);
          magnitude = (G * bodies[i].mass * bodies[j].mass) / (distance * distance);

          factor = magnitude / distance;
          direction.x = bodies[j].position.x - bodies[i].position.x;
          direction.y = bodies[j].position.y - bodies[i].position.y;

          // +force for body i
      #pragma omp atomic
          bodies[i].force.x += factor * direction.x;
      #pragma omp atomic
          bodies[i].force.y += factor * direction.y;

          // -force for body j
      #pragma omp atomic
          bodies[j].force.x -= factor * direction.x;
      #pragma omp atomic
          bodies[j].force.y -= factor * direction.y;
        }
    }
    }

そして、ここでそれを使用する必要がないのはなぜですか：

#pragma omp parallel
{
  vector_t delta_v, delta_p;
  int i;

#pragma omp for
  for (i = 0; i < n_body; i++)
    {
      // calculate delta_v
      delta_v.x = bodies[i].force.x / bodies[i].mass * dt;
      delta_v.y = bodies[i].force.y / bodies[i].mass * dt;

      // calculate delta_p
      delta_p.x = (bodies[i].velocity.x + delta_v.x / 2.0) * dt;
      delta_p.y = (bodies[i].velocity.y + delta_v.y / 2.0) * dt;

       // update body velocity and position
      bodies[i].velocity.x += delta_v.x;
      bodies[i].velocity.y += delta_v.y;
      bodies[i].position.x += delta_p.x;
      bodies[i].position.y += delta_p.y;

      // reset forces
      bodies[i].force.x = bodies[i].force.y = 0.0;


      if (bounce)
    {
      // bounce on boundaries (i.e. it's more like billard)
      if ((bodies[i].position.x < -body_distance_factor) || (bodies[i].position.x > body_distance_factor))
        bodies[i].velocity.x = -bodies[i].velocity.x;
      if ((bodies[i].position.y < -body_distance_factor) || (bodies[i].position.y > body_distance_factor))
        bodies[i].velocity.y = -bodies[i].velocity.y;
    }
    }
}

コードは現在でも機能していますが、その理由がわかりません。手伝って頂けますか？

敬具マイケル

4

2 に答える 2

0

複数のスレッドが同じメモリ位置に書き込む場合、アトミック オペレーターまたはクリティカル セクションを使用して競合状態を防ぐ必要があります。アトミック演算子は高速ですが、より多くの制限があります (たとえば、いくつかの基本的な演算子を使用した POD でのみ動作します) が、あなたの場合はそれらを使用できます。

そのため、スレッドがいつ同じメモリ位置に書き込んでいるかを自問する必要があります。最初のケースでは、外側のループのみを並列化し、内側のループを並列化しiないため、実際には 上のものだけにjアトミック演算子は必要ありません。 ij

最初のケースから例を考えてみましょう。n_body=1014 つのスレッドがあるとします。

Thread one   i =  0-24, j =  1-100, j range = 100
Thread two   i = 25-49, j = 26-100, j range = 75
Thread three i = 50-74, j = 51-100, j range = 50
Thread four  i = 75-99, j = 76-100, j range = 25

まず、各スレッドが同じメモリ位置に書き込みを行っていることがわかります。たとえば、すべてのスレッドは、j=76-100. そのため、アトミック演算子 for が必要ですj。ただし、 と同じメモリ位置に書き込むスレッドはありませんi。そのため、 のアトミック オペレーターは必要ありませんi。

2番目のケースでは、ループが1つしかなく、並列化されているため、同じメモリ位置にスレッドが書き込まれないため、アトミック演算子は必要ありません。

それはあなたの質問に答えますが、コードのパフォーマンスを向上させるための追加のコメントを次に示します。

アトミック演算子とは別の重要な観察事項があります。jスレッド 1 は100 回以上実行されるのに対し、スレッド 4 はj25 回以上実行されることがわかります。したがって、通常はデフォルトのスケジューラである schedule(static) を使用すると、負荷が十分に分散されません。値が大きいほど、n_bodyこれは悪化します。

1 つの解決策は、 schedule(guided)を試すことです。これは以前に使用したことがありませんが、適切なソリューションOpenMP: for scheduleだと思います。「特別な種類の動的スケジューリングがガイドされており、作業が進むにつれて、各タスクにますます小さな反復ブロックが与えられます。」ガイド付きの標準に従って、連続する各ブロックは「number_of_iterations_remaining / number_of_threads」を取得します。だから私たちの例から

Thread one   i =  0-24, j =  1-100, j range = 100
Thread two   i = 25-44, j = 26-100, j range = 75
Thread three i = 45-69, j = 46-100, j range = 55
Thread four  i = 60-69, j = 61-100, j range = 40
Thread one   i = 70-76, j = 71-100
...

スレッドがより均等に分散されていることに注意してください。静的スケジューリングでは、4 番目のスレッドは j で 25 回しか実行されませんでしたが、4 番目のスレッドはj40 回以上実行されます。

しかし、アルゴリズムをもっと注意深く見てみましょう。最初のケースでは、各ボディの重力を計算しています。これを行うには、次の 2 つの方法があります。

//method1 
#pragma omp parallel for
for(int i=0; i<n_body; i++) {
    vector_t force  = 0;
    for(int j=0; j<n_body; j++) {
        force += gravity_force(i,j);
    }
    bodies[i].force = force;
}

しかし、関数gravity_force(i,j)=gravity_force(j,i)なので、2 回計算する必要はありません。だから、より速い解決策を見つけました：

//method2 (see your code and mine below on how to parallelize this)
for(int i=0; i<(n_body-1); i++) {
    for(int j=i+1; j<nbody; j++) {
        bodies[i].force += gravity_force(i,j);
        bodies[j].force += gravity_force(i,j);
    }
}

最初の方法はn_bodies*n_bodies反復を行い、2 番目の方法は (n_bodies-1) nbodies/2 回の反復を行います。最初の順序は n_bodies n_bodies/2 です。 ただし、2 番目のケースは効率的に並列化するのがはるかに困難です (以下のコードと私のコードを参照してください)。Atmoic 操作を使用する必要があり、負荷がバランスしていません。最初の方法は 2 倍の反復を行いますが、負荷は均等に分散され、アトミック操作は必要ありません。両方の方法をテストして、どちらが高速かを確認する必要があります。

2番目の方法を並列化するには、あなたがしたことを行うことができます:

#pragma omp parallel for schedule(static) // you should try schedule(guided)
for(int i=0; i<(n_body-1); i++) {
    for(int j=i+1; j<nbody; j++) {
        //#pragma omp atomic //not necessary on i
        bodies[i].force += gravity_force(i,j);
        #pragma omp atomic   //but necessary on j
        bodies[j].force += gravity_force(i,j);
    }
}

または、より良い解決策は、次のような強制のプライベート コピーを使用することです。

#pragma omp parallel 
{
    vector_t *force = new vector_t[n_body];
    #pragma omp for schedule(static) 
    for (int i = 0; i < n_body; i++) force[i] = 0;
    #pragma omp for schedule(guided)
    for(int i=0; i<(n_body-1); i++) {
        for(int j=i+1; j<nbody; j++) {
            force[i] += gravity_force(i,j);
            force[j] += gravity_force(i,j);
        }
    }
    #pragma omp for schedule(static)
    {
         #pragma omp atomic
         bodies[i].force.x += force[i].x;
         #pragma omp atomic
         bodies[i].force.y += force[i].y;
    }
    delete[] force;
}
于 2013-10-31T14:48:38.903 に答える