OpenMP并行化低效排查及二维Poisson求解代码优化咨询

问题背景 #

尝试用对角线法并行化求解200x200网格的二维Poisson方程，串行版本耗时1分12秒，但OpenMP并行版本无限运行；20x20小网格仅在schedule(dynamic,1024)（近似串行）时可运行。相关代码如下：

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

int xf = -1;
int xl = 1;
int yf = -1;
int yl = 1;

double delta = 0.01;

double q(int i, int j) {
  return (4 - 2 * (pow(xf + j * delta, 2) + pow(yl - i * delta, 2))) *
         pow(delta, 2);
}

double actual(double x, double y, int n) {
  return (pow(x, 2) - 1) * (pow(y, 2) - 1);
}

double error(double **mk, double **new, int n) {
  double sum1 = 0;
  double sum2 = 0;
  double mf;
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      mf = mk[i][j];
      sum1 += pow(mf, 2);
      sum2 += pow(mf - new[i][j], 2);
    }
  }
  return pow((double)sum2 / sum1, 0.5);
}

int main(void) {

  int n = (xl - xf) / delta + 1;

  double **phiactual = (double **)malloc(n * sizeof(double *));
  for (int i = 0; i < n; i++)
    phiactual[i] = (double *)malloc(n * sizeof(double));

  double **phi = (double **)malloc(n * sizeof(double *));
  for (int i = 0; i < n; i++)
    phi[i] = (double *)malloc(n * sizeof(double));

  int i, j, d;

  int iter;
#pragma omp parallel shared(phi, phiactual, delta)
  {
    iter = 0;
#pragma omp for collapse(2)
    for (int q = 0; q < n; ++q) {
      for (int w = 0; w < n; ++w) {
        phiactual[q][w] = actual(xf + w * delta, yl - q * delta, delta);
      }
    }

    while (error(phiactual, phi, n) > 0.01) {
      for (int g = 0; g < 2 * n - 5; ++g) {
        if (g < n - 3) {
          i = 1;
          j = g + 1;
          d = g + 1;
        } else {
          i = g - n + 4;
          j = n - 2;
          d = 2 * n - 5 - g;
        }

#pragma omp for
        for (int k = 0; k < d; ++k) {
          phi[i][j] = 0.25 * (((double)(phi[i + 1][j] + phi[i][j + 1] +
                                        phi[i - 1][j] + phi[i][j - 1])) +
                              q(i, j));
          ++i;
          --j;
        }
      }
      iter++;
    }
#pragma omp single
    printf("%i\n", iter);
  }

  for (int i = 0; i < n; i++)
    free(phi[i]);
  free(phi);

  for (int i = 0; i < n; i++)
    free(phiactual[i]);
  free(phiactual);

  return 0;
}

并行化无限运行/低效的核心原因 #

• 数据竞争破坏迭代依赖：对角线法要求每个点的更新基于上一轮迭代的相邻点值，但当前代码直接在共享的phi数组上原地更新，多个线程同时读写同一内存位置，导致数据竞争。迭代过程中，后续点可能读取到线程刚写入的当前轮次值，而非上一轮的正确值，完全破坏收敛逻辑，导致迭代无法终止。小网格用大动态块时，线程执行区域几乎不重叠，冲突极少，因此能勉强运行，但失去并行加速效果。
• 误差计算的线程安全问题：error函数在并行区域的while循环中被调用，函数内的sum1、sum2是局部变量，但多个线程同时遍历共享的phiactual和phi数组时，没有任何同步机制，会导致误差计算结果错误，error(...) > 0.01的判断永远成立，陷入无限循环。
• 并行区域变量管理错误：iter变量在并行区域内被所有线程同时初始化iter=0，虽然最后用#pragma omp single输出，但迭代计数的更新没有同步，属于未定义行为；同时i、j、d是并行区域内的共享变量，在对角线遍历的循环中被多个线程修改，进一步加剧数据竞争。

串行代码效率优化建议 #

• 改用连续内存存储二维数组：当前二级指针的动态分配会导致内存碎片化，缓存命中率极低。改为单块连续内存的一维数组，例如：

  double *phi = (double*)malloc(n*n*sizeof(double));
  // 访问方式：phi[i*n + j]
  

  连续内存能大幅提升缓存利用率，减少内存访问开销。
• 替换高开销数学函数：pow(x,2)的开销远大于直接计算x*x，将q、actual、error函数中的所有平方运算替换为直接乘法，避免不必要的函数调用。
• 预计算网格坐标：提前计算所有网格点的x、y坐标并存储到数组中，避免每次调用q和actual时重复计算xf + j*delta、yl - i*delta这类表达式，减少计算量。
• 初始化迭代数组：当前phi数组未初始化，初始值为随机垃圾值，导致迭代次数大幅增加。应将phi初始化为边界条件或0，减少迭代步数。
• 优化误差计算逻辑：将error函数中的pow(mf - new[i][j], 2)替换为直接平方运算；若允许，可改用最大绝对误差替代相对误差，减少求和运算的开销，加快收敛判断速度。

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