如何对删除重叠球体的代码进行OpenMP并行化处理？

• 变量私有逻辑错误：你将只读的坐标vectorxcentro/ycentro/zcentro声明为private，OpenMP会为每个线程创建独立的空vector副本，访问空容器的下标直接触发段错误，这类只读全局容器应设为shared，或者不用显式声明，并行块外定义的变量默认是shared属性。
• 嵌套并行指令错误：两层循环都加#pragma omp parallel for属于错误用法，嵌套并行默认关闭，强制开启会导致线程数超额、调度开销激增，仅需要对外层i循环添加并行指令即可。
• 数据竞争问题：多个线程同时修改r数组的同一个元素时会出现竞争，虽然你此处仅做置0操作，多次写0最终结果一致，但如果是其他写入操作会直接导致结果错误；另外你没有跳过已标记为0的球体，大量无效计算拉低效率。
• 边界逻辑错误：循环终止条件写为i < r.size()-1、j < r.size()-1，会漏掉最后一个球体的判断。

首先你可以先跳过r[i]已经为0的外层循环，减少无效计算，OpenMP仅需要对外层循环做并行，私有内层循环变量j和距离变量即可，修改后代码如下：

vector<double> xcentro, ycentro, zcentro, r;
r.reserve(34000000);
xcentro.reserve(34000000);
ycentro.reserve(34000000);
zcentro.reserve(34000000);
// ... 读取文件填充数据 ...

#pragma omp parallel for private(j, dist_square)
for (size_t i = 0; i < r.size(); i++)
{
    if (r[i] == 0) continue; // 已经失效的球体直接跳过
    for (size_t j = i + 1; j < r.size(); j++)
    {
        if (r[j] == 0) continue; // 已经失效的球体直接跳过
        double dx = xcentro[i] - xcentro[j];
        double dy = ycentro[i] - ycentro[j];
        double dz = zcentro[i] - zcentro[j];
        auto dist_square = dx*dx + dy*dy + dz*dz;
        auto r_sum = r[i] + r[j];
        if (dist_square < r_sum * r_sum)
        {
            r[j] = 0;
            r[i] = 0;
        }
    }
}

注：此处仅做置0操作，多线程同时对同一个元素写0不会出现结果异常，无需加锁，不会影响最终结果。

你当前的O(n²)算法仅适合万级以下的球体计算，百万级球体的计算量会达到1e12次，就算并行也无法在合理时间内跑完，建议引入空间划分优化：

• 将三维空间划分为边长等于最大球体直径的网格，每个球体仅存入它所在的网格中
• 每个球体仅需要和同网格、相邻26个网格内的球体做重叠判断，可将计算量降低到O(n)级别

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