OpenMP并行Jacobi热传导程序速度慢于串行代码问题排查

项目背景 #

我是OpenMP初学者，需要完成基于OpenMP的热传导问题求解项目：采用Jacobi迭代法求解二维矩阵，模拟四周边壁温度固定的平板内部未知温度分布。给定基础代码后需要完成三项任务：

• 统计串行代码运行时长
• 并行化串行代码并开展性能对比
• 尽可能对并行代码做进一步优化

遇到的问题 #

运行串行代码并完成初步并行改造后开展对比测试，发现未做针对性优化时，串行代码性能始终优于并行版本；分别测试100×100、300×300等矩阵规模均得到相同结果，且线程数越高程序运行速度越慢。已知小规模问题下并行线程开销会抵消并行收益，但所选测试规模已足够大，不确定当前并行实现是否存在明显错误。

运行环境 #

8核Apple Mac M1设备，已注释非必要打印语句排除输出开销影响，分别测试串行版本、2线程、4线程、6线程并行版本，并行版本性能均不及串行版本。

附：代码与测试结果 #

串行代码 #

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

int main(int argc, char *argv[])
{

    int m; 
    int n;
    double tol;// = 0.0001;

    int i, j, iter;

    m = atoi(argv[1]);
    n = atoi(argv[2]);
    tol = atof(argv[3]);

    /**
     * @var t, tnew, diff, diffmax,
     * t is the old temprature array,tnew is the new array
     */
    double t[m+2][n+2], tnew[m+1][n+1], diff, difmax;

    /**
     * Timer variables
     * @var start, end 
     */
    double start, end;

    printf("%d %d %lf\n",m,n, tol);

    start = omp_get_wtime();

    // initialise temperature array
    for (i=0; i <= m+1; i++) {
        for (j=0; j <= n+1; j++) {
            t[i][j] = 30.0;
        }
    }

    // fix boundary conditions
    for (i=1; i <= m; i++) {
        t[i][0] = 33.0;
        t[i][n+1] = 42.0;
    }
    for (j=1; j <= n; j++) {
        t[0][j] = 20.0;
        t[m+1][j] = 25.0;
    }

    // main loop
    iter = 0;
    difmax = 1000000.0;
    while (difmax > tol) {

        iter++;

        // update temperature for next iteration
        for (i=1; i <= m; i++) {
            for (j=1; j <= n; j++) {
                tnew[i][j] = (t[i-1][j] + t[i+1][j] + t[i][j-1] + t[i][j+1]) / 4.0;
            }
        }

        // work out maximum difference between old and new temperatures
        difmax = 0.0;
        for (i=1; i <= m; i++) {
            for (j=1; j <= n; j++) {
                diff = fabs(tnew[i][j]-t[i][j]);
                if (diff > difmax) {
                    difmax = diff;
                }
                // copy new to old temperatures
                t[i][j] = tnew[i][j];
            }
        }

    }

    end = omp_get_wtime();

    // print results,
    //Loop tempratures commented out to save performance
    printf("iter = %d  difmax = %9.11lf\n", iter, difmax);
    printf("Time in seconds: %lf \n", end - start);
}

初步并行代码 #

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

int main(int argc, char *argv[])
{
    int m; 
    int n;
    double tol;// = 0.0001;

    int i, j, iter, numThreads;

    m = atoi(argv[1]);
    n = atoi(argv[2]);
    tol = atof(argv[3]);
    numThreads = atoi(argv[4]);

    double t[m+2][n+2], tnew[m+1][n+1], diff, diffmax, privDiffmax;

    double start, end;

    omp_set_num_threads(numThreads);
    start = omp_get_wtime();

    // 初始化并行区域
    #pragma omp parallel private(i, j)
    {
        #pragma omp for collapse(2) schedule(static)
        for (i=0; i <= m+1; i++) {
            for (j=0; j <= n+1; j++) {
                t[i][j] = 30.0;
            }
        }

        // 边界条件赋值
        #pragma omp for schedule(static)
        for (i=1; i <= m; i++) {
            t[i][0] = 33.0;
            t[i][n+1] = 42.0;
        }

        #pragma omp for schedule(static)
        for (j=1; j <= n; j++) {
            t[0][j] = 20.0;
            t[m+1][j] = 25.0;
        }

    }   

    // 主迭代循环
    iter = 0;
    diffmax = 1000000.0;

    while (diffmax > tol) {

        iter = iter + 1;

        // 温度更新
        #pragma omp parallel for private(i, j) collapse(2) schedule(static)
        for (i=1; i <= m; i++) {
            for (j=1; j <= n; j++) {
                tnew[i][j] = (t[i-1][j] + t[i+1][j] + t[i][j-1] + t[i][j+1]) / 4.0;
            }
        }

        // 计算最大差值并更新旧数组
        diffmax = 0.0;
        #pragma omp parallel private(i, j, privDiffmax, diff)
        {
            privDiffmax = 0.0;
            #pragma omp for collapse(2) schedule(static)
            for (i=1; i <= m; i++) {
                for (j=1; j <= n; j++) {
                    diff = fabs(tnew[i][j]-t[i][j]);
                    if (diff > privDiffmax) {
                        privDiffmax = diff;
                    }
                    t[i][j] = tnew[i][j];
                }
            }
            #pragma omp critical
            if (privDiffmax > diffmax)
            {
                diffmax = privDiffmax;
            }
        }
    }

    end = omp_get_wtime();

    printf("iter = %d  diffmax = %9.11lf\n", iter, diffmax);
    printf("Time in seconds: %lf \n", end - start);
}

基准测试结果 #

核心问题排查方向 #

并行版本性能反降的核心原因按影响权重排序如下：

1. 迭代内重复启停线程开销过高
   当前代码在while主循环内部每次都使用独立的#pragma omp parallel和#pragma omp parallel for块，意味着每一轮Jacobi迭代都要触发线程组的创建、同步、销毁。Jacobi迭代收敛通常需要上千甚至上万轮，这部分开销会完全盖过并行计算的收益。
   正确做法是把整个while循环放在同一个#pragma omp parallel区域内部，内部循环仅用#pragma omp for做任务划分，避免每轮迭代的线程fork/join开销。

2. collapse(2)破坏缓存局部性
   二维Jacobi迭代本身是按行连续访问的内存模式，直接对外层i循环做并行划分，每个线程负责连续的若干行，内存访问连续，缓存命中率极高；使用collapse(2)把两层循环拍平后，线程的访问区间会被切得很碎，缓存局部性被完全破坏，访问效率远低于串行的顺序访问。
   建议去掉collapse(2)，仅对外层i循环做静态调度，每个线程处理连续行块，最大化L1/L2缓存利用率。

3. 伪共享与栈数组问题
   代码使用变长数组(VLA)在栈上分配温度矩阵，M1架构缓存行大小为64字节，多线程同时写相邻内存地址时会频繁触发缓存行失效，性能骤降。此外定义的tnew数组比t数组小1个维度，存在潜在内存越界风险。
   建议改用malloc在堆上分配数组，同时按缓存行大小做内存对齐，减少跨核缓存同步开销。

另外M1是性能核+能效核的混合架构，6线程调度时任务可能被分配到能效核上，也会导致性能下降，建议优化完上述问题后先测试2、4线程（仅占用性能核）的场景，再逐步扩展线程数。

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