如何用OpenMP并行化C语言贪心算法实现的旅行商问题代码？

一、适合并行化的核心代码段 #

原代码中存在三个高度独立、无数据依赖的计算模块，是并行化的重点：

1. 随机坐标生成：
   生成X、Y数组的循环，每个元素的计算完全独立，无共享数据依赖。

2. 距离矩阵计算：
   双重循环中，每个(i,j)对的距离计算互不干扰（j > i时仅计算一次，然后赋值给对称位置），适合并行处理。

3. 外层起点遍历循环：
   最核心的性能瓶颈——for (int first = 0; first < nn; first++)，每个迭代对应从不同起点出发的贪心路径求解，所有迭代之间无数据依赖（各起点的路径计算完全独立），并行化收益最高。

二、推荐使用的OpenMP构造 #

针对不同模块，选择对应的OpenMP指令：

1. 并行循环（#pragma omp parallel for） #

适用于随机坐标生成、距离矩阵计算、外层起点循环这类无依赖的迭代任务。

2. 线程私有变量（private/firstprivate） #

确保每个线程拥有独立的局部变量副本，避免数据竞争，比如path数组、dist、current等。

3. 临界区（#pragma omp critical） #

用于全局最优解（best和good数组）的更新，保证同一时间只有一个线程修改共享变量，避免数据竞争。

4. 线程局部最优缓存 #

为每个线程维护局部的local_best和local_good，减少全局同步的开销，仅在线程完成所有迭代后合并全局最优。

三、需要规避的潜在陷阱 #

1. 共享变量数据竞争：
   直接在并行循环中修改全局变量best和good会导致结果错误，必须用临界区保护更新操作；或者用线程局部变量缓存最优解，最后一次性合并。

2. 伪随机数线程安全问题：
   原生rand()使用全局状态，多线程调用会导致随机数重复或性能下降，需改用线程安全的rand_r()，每个线程维护独立的随机种子。

3. 提前终止的效率与准确性平衡：
   原代码中if (dist >= best) break;的提前终止逻辑，若直接使用全局best会频繁触发缓存同步，降低性能。建议每个线程维护local_best副本，仅在合并阶段更新全局值，平衡性能和终止时机的准确性。

4. 内存开销问题：
   当nn较大（如5000）时，每个线程需要独立的path和local_good数组，总内存开销为线程数 × nn × 4字节，需确保系统内存足够。

5. 循环负载均衡：
   静态调度（默认）适合nn远大于线程数的场景；若nn较小，可使用schedule(dynamic)动态分配迭代任务，但会增加调度开销，需根据实际情况权衡。

四、优化后的并行代码示例 #

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>
#include <omp.h>

#define N 5000
#define LARGE (2.0f * (N * 10.0f) * (N * 10.0f)) // 用浮点常量替代pow，提升性能

int X[N];
int Y[N];
float distance[N][N + 1];

int main(int na, char * arg[]) {
    assert(na == 2);

    printf("Dimension %s\n", arg[1]);
    int nn = atoi(arg[1]);
    assert(nn <= N);

    // 并行生成随机坐标，使用线程安全的rand_r
    #pragma omp parallel default(none) shared(X, Y, nn) private(i, seed)
    {
        seed = omp_get_thread_num();
        #pragma omp for
        for (int i = 0; i < nn; i++) {
            X[i] = rand_r(&seed) % (nn * 10);
            Y[i] = rand_r(&seed) % (nn * 10);
        }
    }

    // 并行计算距离矩阵
    #pragma omp parallel for default(none) shared(X, Y, distance, nn) private(i, j, dx, dy)
    for (int i = 0; i < nn; i++) {
        distance[i][i] = 0.0f;
        for (int j = i + 1; j < nn; j++) {
            dx = (float)(X[i] - X[j]);
            dy = (float)(Y[i] - Y[j]);
            distance[i][j] = distance[j][i] = sqrtf(dx*dx + dy*dy); // 用sqrtf替代sqrt，提升浮点性能
        }
    }

    float best = LARGE;
    int good[nn];

    // 并行处理外层起点循环，使用线程局部最优缓存
    #pragma omp parallel default(none) shared(best, good, distance, nn) \
        private(first, dist, path, current, i, j, dmin, index) \
        firstprivate(best)
    {
        float local_best = best;
        int local_good[nn];

        #pragma omp for schedule(static) // 静态调度适合大nn场景
        for (int first = 0; first < nn; first++) {
            float dist = 0.0f;
            int path[nn];
            for (int i = 0; i < nn; i++)
                path[i] = -1;
            path[first] = 0;
            int current = first;

            for (int i = 1; i < nn; i++) {
                float dmin = LARGE;
                int index = 0;
                for (int j = 0; j < nn; j++) {
                    if (path[j] == -1 && current != j && distance[current][j] < dmin) {
                        dmin = distance[current][j];
                        index = j;
                    }
                }
                current = index;
                path[current] = i;
                dist += dmin;

                // 使用局部best判断提前终止，减少同步开销
                if (dist >= local_best) {
                    dist = 0.0f;
                    break;
                }
            }

            if (dist > 0.0f) {
                float dmin = distance[current][first];
                dist += dmin;

                // 更新线程局部最优
                if (dist < local_best) {
                    for (int i = 0; i < nn; i++)
                        local_good[path[i]] = i;
                    local_best = dist;
                }

                distance[first][nn] = dist;
            }
        }

        // 临界区合并全局最优解
        #pragma omp critical
        {
            if (local_best < best) {
                best = local_best;
                for (int i = 0; i < nn; i++)
                    good[i] = local_good[i];
            }
        }
    }

    printf("Solution :\n");
    for (int i = 0; i < nn; i++)
        printf("%d\n", good[i]);
    printf("Distance %g == %g\n", best, distance[good[0]][nn]);

    exit(0);
}

五、额外性能优化建议 #

• 用sqrtf替代sqrt、powf替代pow，减少浮点运算的开销（针对单精度浮点）。
• 预计算LARGE为常量值，避免运行时调用pow函数。
• 对于超大nn，可以考虑使用动态调度schedule(dynamic, chunk_size)，根据硬件调整chunk_size平衡负载。

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