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

3

两个代码示例中的第二个，循环的每个并行迭代都在数组的元素 [i] 上工作，并且从不查看任何相邻元素。因此，循环的每次迭代对循环的任何其他迭代都没有影响，并且它们都可以同时执行而无需担心。

然而，在第一个代码示例中，循环的每个并行迭代都可以使用索引 [j] 读取和写入 body 数组中的任何位置。这意味着两个线程可能试图同时更新相同的内存位置，或者一个线程可能正在写入另一个线程正在读取的位置。为了避免竞争条件，您需要确保写入是原子的。

于 2013-09-16T18:13:10.703 回答
0

当多个线程写入同一内​​存位置时，您需要使用原子运算符或临界区来防止竞争条件。原子运算符速度更快但有更多限制（例如，它们仅在具有一些基本运算符的 POD 上运行）但在您的情况下您可以使用它们。

因此，您必须问自己何时线程正在写入相同的内存位置。在第一种情况下，您只并行化外循环i而不是内循环，j因此您实际上并不需要i仅在j.

让我们考虑第一个案例中的一个例子。假设n_body=101有 4 个线程。

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.

在第二种情况下，您只有一个循环并且它是并行化的，因此没有线程写入相同的内存位置，因此您不需要原子运算符。

这回答了你的问题，但这里有一些额外的评论可以提高你的代码的性能：

还有另一个独立于原子运算符的重要观察。您可以看到线程 1 运行超过j100 次，而线程 4 仅运行超过j25 次。因此，使用通常是默认调度程序的 schedule(static) 不能很好地分配负载。对于较大的值，n_body这将变得更糟。

一种解决方案是尝试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
...

现在请注意，线程分布得更均匀。使用静态调度，第四个线程只运行了 j 25 次，而现在第四个线程运行了j40 次。

但是让我们更仔细地看看你的算法。在第一种情况下，您正在计算每个物体上的重力。这里有两种方法可以做到这一点：

//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)所以你不需要计算两次。所以你找到了一个更快的解决方案：

//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进行迭代，第二种方法进行 (n_bodies-1) nbodies/2 次迭代，一阶是 n_bodies n_bodies/2。 但是，第二种情况更难以有效地并行化（请参阅下面的代码和我的代码）。它必须使用 atmoic 操作并且负载不平衡。第一种方法执行两倍的迭代，但负载分布均匀，不需要任何原子操作。您应该测试这两种方法，看看哪一种更快。

要并行化第二种方法，您可以执行以下操作：

#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 回答