Увеличение производительности программы n-body с использованием OpenMP

Question

Моя цель - повысить производительность кода, имитирующего проблему n-тела.

Здесь рассчитывается время.Двумя функциями, которые необходимо распараллелить, являются функции convert_forces () и * move_body () , но так как переменная управления цикла t - это двойное число, у меня не может быть # параллельной прагмы для оператора .

t0 = gettime ();
  for (t = 0; t < t_end; t += dt)
    {
      // draw bodies
      show_bodies (window);

      // computation
      calculate_forces ();
      move_bodies ();
    }

  // print out calculation speed every second
  t0 = gettime () - t0;

Две функции execute_forces () и move_body () с соответствующими директивами, которые я использовал, следующие:

static void
calculate_forces ()
{
  double distance, magnitude, factor, r;
  vector_t direction;
  int i, j;
 #pragma omp parallel private(distance,magnitude,factor,direction)
 {
   #pragma omp for private(i,j)
   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 critical
         {
          bodies[i].force.x += factor * direction.x;
          bodies[i].force.y += factor * direction.y;

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

}
static void
move_bodies ()
{
  vector_t delta_v, delta_p;
  int i;
 #pragma omp parallel private(delta_v,delta_p,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
      #pragma omp critical
      {
      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;
        }
    }
}

ЗначенияТел. скорость и тел.позиция меняются в функции перемещения тел, но я не смог использовать сокращение.

Существует также функция контрольной суммы для вычисления, равна ли вычисленная контрольная сумма контрольной сумме.Эта функция выглядит следующим образом:

static unsigned long
checksum()
{
  unsigned long checksum = 0;

  // initialize bodies
  for (int i = 0; i < n_body; i++)
    {
      // random position vector
      checksum += (unsigned long)round(bodies[i].position.x);
      checksum += (unsigned long)round(bodies[i].position.y);
    }

  return checksum;
}

Эта функция использует ранее вычисленные значения body.position.x и body.position.y, которые были вычислены в функции move_bodies, поэтому я и использовала критический блок.при расчете тех значений, которые не дают правильного ответа.Кто-нибудь может дать мне некоторое представление о том, где я иду не так?Заранее спасибо.