OpenMP 锁与关键锁

Question

我正在尝试使用锁和关键部分来确保循环线程安全。这是代码：

#pragma omp parallel for num_threads(4) private(k, f_part_k, len, len_3, mg, fact)
for (k = part+1; k < n; k++) {
  /* Compute force on part due to k */
  f_part_k[X] = curr[part].s[X] - curr[k].s[X];
  f_part_k[Y] = curr[part].s[Y] - curr[k].s[Y];
  len = sqrt(f_part_k[X]*f_part_k[X] + f_part_k[Y]*f_part_k[Y]);
  len_3 = len*len*len;
  mg = -G*curr[part].m*curr[k].m;
  fact = mg/len_3;
  f_part_k[X] *= fact;
  f_part_k[Y] *= fact;

  /* Add force in to total forces */
  omp_set_lock(&(locks[k]));
  //#pragma omp critical
  {
      forces[part][X] += f_part_k[X];
      forces[part][Y] += f_part_k[Y];
      forces[k][X] -= f_part_k[X];
      forces[k][Y] -= f_part_k[Y];
  }
  omp_unset_lock(&(locks[k]));
}

for (i = 0; i < n; i++)
    omp_destroy_lock(&(locks[i]));
} 

当我仅使用被注释掉的关键指令时，结果很好，即与顺序版本的结果相匹配。但是，如果我使用代码中所示的锁，结果会很糟糕。我想我误解了锁的概念，因为根据我的理解，使用这种锁方法对forces数组的写访问应该是安全的。你能指出我正确的方向吗？

Answer 1

我认为您的代码的问题是以下方面的竞争条件：

omp_set_lock(&(locks[k]));
{
    forces[part][X] += f_part_k[X]; // Race condition for different k
    forces[part][Y] += f_part_k[Y]; // Race condition for different k
    forces[k][X] -= f_part_k[X]; 
    forces[k][Y] -= f_part_k[Y]; 
}
omp_unset_lock(&(locks[k]));

事实上，对于 的不同值k，多个线程尝试写入forces[part][X]和forces[part][Y]。此外，我认为不需要显式同步对forces[k][X]和 的访问forces[k][Y]，因为每个线程都会更新自己的k。

如果您想尝试提供正确语义的不同同步结构，您可以尝试：

原子级同步

#pragma omp atomic
forces[part][X] += f_part_k[X];
#pragma omp atomic
forces[part][Y] += f_part_k[Y];

forces[k][X] -= f_part_k[X]; 
forces[k][Y] -= f_part_k[Y]; 

显式锁定

omp_set_lock(&lock);
{
  forces[part][X] += f_part_k[X];
  forces[part][Y] += f_part_k[Y];
}
omp_unset_lock(&lock);

forces[k][X] -= f_part_k[X]; 
forces[k][Y] -= f_part_k[Y];

命名临界区

#pragma omp critical(PART)
{
  forces[part][X] += f_part_k[X];
  forces[part][Y] += f_part_k[Y];
}
forces[k][X] -= f_part_k[X]; 
forces[k][Y] -= f_part_k[Y]; 

critical我建议您阅读和atomic结构的定义（第 2.8.2 和 2.8.5节），并查看示例A.19.1c、A.22.*和A.45.1c

也就是说，在您提出的情况下，我会尝试以下操作：

float fredx = 0.0f;
float fredy = 0.0f;
#pragma omp parallel for private(k, f_part_k, len, len_3, mg, fact) reduction(+:fredx,fredy)
for (k = part+1; k < n; k++) {
  /* Compute force on part due to k */
  f_part_k[X] = curr[part].s[X] - curr[k].s[X];
  f_part_k[Y] = curr[part].s[Y] - curr[k].s[Y];
  len = sqrt(f_part_k[X]*f_part_k[X] + f_part_k[Y]*f_part_k[Y]);
  len_3 = len*len*len;
  mg = -G*curr[part].m*curr[k].m;
  fact = mg/len_3;
  f_part_k[X] *= fact;
  f_part_k[Y] *= fact;

  /* Add force in to total forces */
  fredx += f_part_k[X];
  fredy += f_part_k[Y];

  forces[k][X] -= f_part_k[X];
  forces[k][Y] -= f_part_k[Y];            
}

forces[part][X] += fredx;
forces[part][Y] += fredy;

以避免任何显式同步。