使用OpenMP并行化C语言嵌套循环遇阻，寻求解决方案

我尝试用OpenMP并行化C语言里的这个函数，但一直没成功。我知道迭代过程和两行代码之间有依赖关系，但卡壳了，希望得到帮助。

串行参考代码 #

void loop_reference(double C[N][N], double A[N][N], double B[N][N]) {
  size_t i, j;

  for (i = 1; i < N; i++) {
    for (j = 1; j < N; j++) {
      B[i][j] = B[i-1][j] + A[i][j-1];
      A[i][j] = A[i-1][j] + C[i][j];
    }
  }
}

我的尝试（无效） #

我试过把两行代码拆分后分别并行化，但没用，尝试的并行代码如下：

void loop_parallel(double C[N][N], double A[N][N], double B[N][N]) {
  size_t i, j;

  #pragma omp parallel for private(i, j) shared(A, B, C)
  for (i = 1; i < N; i++) {
    for (j = 1; j < N; j++) {
      A[i][j] = A[i-1][j] + C[i][j];
      
    }
  }

  #pragma omp parallel for private(i, j) shared(A, B, C)
  for (i = 1; i < N; i++) {
    for (j = 1; j < N; j++) {
      A[i][j] = A[i-1][j] + C[i][j];
    }
  }
}

问题分析 #

你的并行代码的错误 #

1. 完全遗漏了B数组的计算逻辑，两次循环都只重复计算A，直接导致结果和串行代码不一致
2. 强行拆分计算顺序，原代码是每个(i,j)先更新B再更新A，而你的代码先批量更新所有A，彻底破坏了原有的依赖关系

原代码的依赖关系 #

要并行化必须先理清依赖：

• 行内j维度依赖：同一行i中，B[i][j]依赖A[i][j-1]（同行前一列的A），而A[i][j-1]是同一行j-1位置刚计算的结果，因此j循环必须顺序执行，无法直接并行
• 行间i维度依赖：B[i][j]依赖B[i-1][j]（上一行同列的B），A[i][j]依赖A[i-1][j]（上一行同列的A），因此i循环也无法直接并行

可行的并行方案：波前并行 #

虽然i和j循环无法直接并行，但可以采用波前并行（按i+j=常数的对角线分组），同一对角线上的所有(i,j)之间没有依赖关系，可以安全并行。

并行代码示例 #

void loop_parallel(double C[N][N], double A[N][N], double B[N][N]) {
  size_t i, j, k;
  // 波前的对角线k = i + j，范围从2（i=1,j=1）到2*(N-1)（i=N-1,j=N-1）
  #pragma omp parallel for private(i, j) shared(A, B, C) schedule(dynamic)
  for (k = 2; k <= 2*(N-1); k++) {
    // 遍历当前对角线上的所有(i,j)
    for (i = 1; i < N; i++) {
      j = k - i;
      if (j >= 1 && j < N) {
        B[i][j] = B[i-1][j] + A[i][j-1];
        A[i][j] = A[i-1][j] + C[i][j];
      }
    }
  }
}

说明 #

• 外层循环按对角线k=i+j遍历，同一k的所有(i,j)计算不需要彼此的结果，因此可以并行
• 使用schedule(dynamic)可以让负载更均衡，避免不同对角线的元素数量差异导致的负载不均
• 变量i,j设为private，A,B,C设为shared，符合OpenMP的变量作用域要求

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