OpenMP并行化含Eigen类成员的粒子箱更新代码时原子操作报错求助

我正在进行N个粒子在指定范围内相互作用的模拟。为避免粒子的N²复杂度计算，我对粒子进行空间排序，将粒子索引存入数组，使每个粒子指向同箱内的另一个粒子。我已编写C++串行代码，现尝试实现OpenMP版本以支持更多粒子。

类定义 #

class Boxes
{
    int m_NX;
    int m_NY;
    int m_Nboxes;
    Eigen::ArrayXi m_boxes;
    // ...
};

class Particles
{
    int m_nbParticles;
    Eigen::ArrayXd m_positions;
    Eigen::ArrayXi m_nextParticles;
    // ...
};

串行排序代码 #

void updateBoxes(Boxes &p_boxes, Particles &p_particles)
{
    // ...
    for (int i = 0; i < p_particles.m_nbParticles; i++)
    {
        int indexX = p_particles.position(i).x() / dX;
        int indexY = p_particles.position(i).y() / dX;
        int indexBox = indexX + NXboxes*indexY;
        p_particles.m_nextParticles[i] = p_boxes.m_boxes[indexBox];
        p_boxes.m_boxes[indexBox] = i;
    }
}

并行化尝试与错误 #

我尝试添加OpenMP指令并行化，但编译报错：

#pragma omp parallel for
for (int i = 0; i < p_particles.size(); i++)
{
    int indexX = p_particles.position(i).x() / dX;
    int indexY = p_particles.position(i).y() / dX;
    int indexBox = indexX + NXboxes*indexY;
    #pragma omp atomic
    p_particles.m_nextParticles[i] = p_boxes.m_boxes[indexBox];
    #pragma omp atomic
    p_boxes.m_boxes[indexBox] = i;
}

错误信息：

error: the statement for 'atomic' must be an expression statement of form '++x;', '--x;', 'x++;', 'x--;', 'x binop= expr;', 'x = x binop expr' or 'x = expr binop x', where x is an l-value expression with scalar type

这段代码单线程下占总耗时约8%，但线程数越多占比越高，请问如何正确并行化？

你用atomic的方式不符合OpenMP规范，因为atomic只支持错误提示中列出的特定操作（自增、自减、复合赋值等），普通赋值不在其支持范围内。这段代码的核心是每个箱子的链表头更新必须是原子的读-改-写操作：先读取当前箱子的头粒子索引，赋值给当前粒子的nextParticles，再将箱子头更新为当前粒子索引，两步必须原子完成，否则多线程操作同一箱子会导致链表断裂。

方案1：使用OpenMP临界区（简单直接） #

给每个箱子的操作加临界区，确保同一箱子的修改不会被并发执行：

void updateBoxes(Boxes &p_boxes, Particles &p_particles)
{
    // ...
    #pragma omp parallel for
    for (int i = 0; i < p_particles.m_nbParticles; i++)
    {
        int indexX = p_particles.position(i).x() / dX;
        int indexY = p_particles.position(i).y() / dX;
        int indexBox = indexX + NXboxes * indexY;

        #pragma omp critical(box_update)
        {
            p_particles.m_nextParticles[i] = p_boxes.m_boxes[indexBox];
            p_boxes.m_boxes[indexBox] = i;
        }
    }
}

• 给临界区命名（如box_update），避免与代码中其他临界区冲突。
• 如果箱子数量远大于线程数，临界区的性能开销会很小；若箱子数量少，建议用方案2。

方案2：使用原子交换操作（性能更优） #

利用原子交换指令一步完成“读取旧值+写入新值”的原子操作，避免临界区的锁开销：

void updateBoxes(Boxes &p_boxes, Particles &p_particles)
{
    // ...
    #pragma omp parallel for
    for (int i = 0; i < p_particles.m_nbParticles; i++)
    {
        int indexX = p_particles.position(i).x() / dX;
        int indexY = p_particles.position(i).y() / dX;
        int indexBox = indexX + NXboxes * indexY;

        // 原子交换：将i写入m_boxes[indexBox]，返回原来的箱子头索引
        int old_head = __atomic_exchange_n(&p_boxes.m_boxes[indexBox], i, __ATOMIC_SEQ_CST);
        p_particles.m_nextParticles[i] = old_head;
    }
}

• __atomic_exchange_n是GCC/Clang的内置原子函数，__ATOMIC_SEQ_CST确保内存顺序的一致性，完全满足粒子模拟的需求。
• 若编译器支持OpenMP 5.0+，也可以用标准OpenMP原子交换语法替代内置函数：

  int old_head;
  #pragma omp atomic exchange
  old_head = p_boxes.m_boxes[indexBox];
  p_boxes.m_boxes[indexBox] = i;
  p_particles.m_nextParticles[i] = old_head;
  

额外优化提示 #

• 确保dX、NXboxes等变量是并行区域内可见的共享变量，必要时用firstprivate声明。
• 粒子位置读取（position(i).x()）是只读操作，Eigen的ArrayXd在只读场景下线程安全，无需额外同步。
• 若箱子数量极少，可考虑先按箱子分组粒子，再串行处理每组箱子，减少并发冲突。

内容的提问来源于stack exchange，提问作者Hunken