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]));
} 

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

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 构造 here 的定义（第 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;

避免任何显式同步。