如何使用OpenMP并行化累积概率函数的计算？

原代码存在的核心问题 #

• 时间复杂度为O(n²)，嵌套循环在数据量大时运行效率极低
• OpenMP使用错误：未实际开启循环并行，额外添加的atomic指令反而增加了运行开销，且atomic write不适用累加场景
• 存在语法错误：重复定义变量aux、vector声明末尾误用逗号、循环终止条件!= r.size()-1会漏掉最后一个元素
• 浮点数直接用==比较，存在精度误差风险

第一阶段：OpenMP 并行化改造（针对原有嵌套循环逻辑） #

原逻辑中每个i对应的计算完全独立，无数据依赖，只需要对外层循环开启并行，明确声明私有变量即可，无需任何atomic操作，改造后可直接获得和CPU核心数近似线性的性能提升：

#include <vector>
#include <cmath>
#include <omp.h>

// 调用前需提前给contador、contador1、vol分配和r相同大小的空间
void calc_cdf(std::vector<double>& r, std::vector<double>& contador, 
              std::vector<double>& contador1, std::vector<double>& vol) {
    int n = r.size();
    // 对外层i循环开启并行，私有变量自动分配给每个线程，无数据竞争
    #pragma omp parallel for default(none) shared(r, contador, contador1, vol, n)
    for (int i = 0; i < n; i++) {
        double aux = r[i];
        double contar = 0, contar1 = 0, vol_i = 0.0;
        const double pi = 3.141592653589793;
        for (int j = 0; j < n; j++) {
            if (aux <= r[j]) {
                contar++;
                vol_i += 4.0 * pi * pow(r[j], 3) / 3.0;
            }
            // 浮点数比较增加精度阈值，避免误差
            if (fabs(aux - r[j]) < 1e-9) {
                contar1++;
            }
        }
        contador[i] = contar;
        contador1[i] = contar1;
        vol[i] = vol_i;
    }
}

第二阶段：算法复杂度优化（性能提升幅度远大于纯并行） #

原嵌套循环逻辑可以通过排序+二分查找+后缀和优化到O(nlogn)，数据量越大提升越明显，优化后也可以配合OpenMP进一步加速：

1. 提前计算所有半径对应的体积，生成半径-体积配对数组
2. 对配对数组按半径升序排序
3. 预计算体积的后缀和数组，后缀和数组第k位表示从第k位到末尾的体积总和
4. 对每个原半径值，用二分查找找到排序后第一个大于等于该值的位置，直接得到符合条件的元素数量和对应体积总和

优化后代码示例：

#include <vector>
#include <algorithm>
#include <cmath>
#include <omp.h>

void calc_cdf_fast(std::vector<double>& r, std::vector<double>& contador, 
                   std::vector<double>& contador1, std::vector<double>& vol) {
    int n = r.size();
    std::vector<std::pair<double, double>> r_vol(n);
    const double pi = 3.141592653589793;
    // 预计算每个半径对应的体积，可并行
    #pragma omp parallel for default(none) shared(r, r_vol, pi, n)
    for (int i = 0; i < n; i++) {
        r_vol[i].first = r[i];
        r_vol[i].second = 4.0 * pi * pow(r[i], 3) / 3.0;
    }
    // 按半径升序排序
    std::sort(r_vol.begin(), r_vol.end());
    // 预计算体积后缀和
    std::vector<double> suffix_sum(n + 1, 0.0);
    for (int i = n - 1; i >= 0; i--) {
        suffix_sum[i] = suffix_sum[i + 1] + r_vol[i].second;
    }
    // 对每个原半径查询结果，可并行
    #pragma omp parallel for default(none) shared(r, r_vol, suffix_sum, contador, contador1, vol, n)
    for (int i = 0; i < n; i++) {
        double aux = r[i];
        // 二分查找边界
        auto low = std::lower_bound(r_vol.begin(), r_vol.end(), std::make_pair(aux, 0.0));
        auto high = std::upper_bound(r_vol.begin(), r_vol.end(), std::make_pair(aux, 1e18));
        contador[i] = n - (low - r_vol.begin());
        contador1[i] = high - low;
        vol[i] = suffix_sum[low - r_vol.begin()];
    }
}

该版本在数据量大于1000时，性能比原嵌套循环高几十到上百倍，配合OpenMP可以进一步获得数倍提升。

编译注意事项 #

• 编译时需要添加OpenMP编译参数：GCC/Clang加-fopenmp，MSVC加/openmp
• 如果要使用标准库的M_PI常量，需要在包含cmath前定义#define _USE_MATH_DEFINES

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