OpenMP并行化基数排序C代码的循环并行性及指令子句咨询

两个待判断循环的并行性分析 #

1. count数组前缀和循环 #

for (i = 1; i < 10; i++)
    count[i] += count[i - 1];

这个循环完全不需要并行，原因有两点：

• 迭代次数仅9次，OpenMP创建线程、同步的开销远大于计算本身的收益，并行反而会拖慢速度
• 迭代存在真数据依赖：第i次迭代的计算依赖第i-1次的count结果，强行并行需要实现额外的并行前缀求和逻辑，对于长度仅为10的数组来说完全得不偿失
  直接保留当前串行写法即可，不会有性能损失。

2. 倒序构建output数组的循环 #

你当前的写法存在严重问题，且并行效率极低：

• 仅给output的赋值加了atomic write，但下方的count[(arr[i] / exp) % 10]--是读-改-写操作，没有原子保护，会出现数据竞争，最终排序结果错误
• 就算给count的递减也加上原子操作，多个线程同时修改count数组的10个元素会产生极高的竞争开销，并行收益会被完全抵消，甚至比串行运行还慢
  这个循环不适合用当前加原子的方式并行，需要调整整体count统计逻辑来避免竞争。

现有代码的其他问题 #

• 统计count的第一个并行循环用atomic递增count元素，同样存在严重的多线程竞争，原子操作开销极大
• 用栈上可变长数组unsigned output[n]，当n较大时会触发栈溢出，建议改用malloc动态分配
• 循环变量i用int类型，当n大于INT_MAX时会溢出，建议统一用size_t类型
• 采用10为基数的排序策略，对于32/64位无符号数来说需要的迭代次数过多，效率低于256基数（8位一个分段）的实现。

优化后的并行实现方案 #

核心思路是用线程局部count数组避免统计阶段的原子竞争：

1. 统计阶段：每个线程维护自己的局部count[10]数组，并行统计自己负责的那部分数据的桶计数，完全无竞争不需要原子操作
2. 合并阶段：串行把所有线程的局部count数组合并为全局count数组，仅需要10 * 线程数次操作，开销极低
3. 前缀和阶段：保留原来的串行前缀和即可
4. 写output阶段：用atomic capture原子获取count递减后的值，保证读写一致性，竞争开销远低于原写法的两次独立原子操作

优化后countSort核心代码示例 #

void countSort(size_t n, unsigned arr[n], unsigned exp) {
    unsigned *output = malloc(n * sizeof(unsigned));
    int global_count[10] = {0};
    int num_threads = omp_get_max_threads();
    // 分配每个线程的局部count数组
    int (*local_count)[10] = malloc(num_threads * sizeof(*local_count));
    memset(local_count, 0, num_threads * sizeof(*local_count));

    // 并行统计局部count，无原子操作无竞争
    #pragma omp parallel default(none) shared(arr, n, exp, local_count)
    {
        int tid = omp_get_thread_num();
        #pragma omp for
        for (size_t i = 0; i < n; i++) {
            int digit = (arr[i] / exp) % 10;
            local_count[tid][digit]++;
        }
    }

    // 合并局部count到全局count
    for (int t = 0; t < num_threads; t++) {
        for (int d = 0; d < 10; d++) {
            global_count[d] += local_count[t][d];
        }
    }

    // 串行前缀和
    for (int i = 1; i < 10; i++) {
        global_count[i] += global_count[i-1];
    }

    // 倒序写output，用atomic capture保证count操作的原子性
    #pragma omp parallel for
    for (size_t i = n-1; i < n; i--) { // size_t是无符号类型，i减到0后再减1会变为SIZE_MAX，自动退出循环
        int digit = (arr[i] / exp) % 10;
        int pos;
        // 原子捕获count递减后的值，一次原子操作完成读改写
        #pragma omp atomic capture
        pos = --global_count[digit];
        output[pos] = arr[i];
    }

    // 拷贝回原数组
    #pragma omp parallel for
    for (size_t i = 0; i < n; i++) {
        arr[i] = output[i];
    }

    free(output);
    free(local_count);
}

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