c++ - 是否可以并行化这个 for 循环？

Question

我得到了一些使用 OpenMP 进行并行化的代码，在各种函数调用中，我注意到这个 for 循环在计算时间上有一些好处。

  double U[n][n];
  double L[n][n];
  double Aprime[n][n];
  for(i=0; i<n; i++) {
    for(j=0; j<n; j++) {
      if (j <= i) {
          double s;
          s=0;
          for(k=0; k<j; k++) {
            s += L[j][k] * U[k][i];
          } 
          U[j][i] = Aprime[j][i] - s;
      } else if (j >= i) {
          double s;
          s=0;
          for(k=0; k<i; k++) {
            s += L[j][k] * U[k][i];
          }
          L[j][i] = (Aprime[j][i] - s) / U[i][i];
      }
    }

然而，在尝试将其并行化并在各处应用一些信号量之后(没有运气)，我开始意识到 else if 条件对早期的 有很强的依赖性if (L[j][i] 是一个用U[i][i] 处理过的数字，可以在早期的 if)，在我看来，由于竞争条件，它是不可并行化的。

是否可以并行化此代码，使 else if 仅在较早的 if 已经完成时执行？

Answer 1

在尝试并行化之前，先尝试简化。

例如，if可以完全消除。

此外，代码访问矩阵的方式会导致最差 缓存性能。这可能是真正的瓶颈。

注意:在下面的更新 #3 中，我做了基准测试和缓存友好版本 fix5 ，从更新 #2 开始，性能比原来的高出 3.9 倍。

我已经分阶段进行了清理，因此您可以看到代码转换。

有了这个，应该可以添加 omp指令成功。正如我在顶部评论中提到的，变量的全局范围与函数范围会影响可能需要的更新类型(例如 omp atomic update 等)

---

作为引用，这是您的原始代码:

double U[n][n];
double L[n][n];
double Aprime[n][n];

for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
        if (j <= i) {
            double s;

            s = 0;
            for (k = 0; k < j; k++) {
                s += L[j][k] * U[k][i];
            }
            U[j][i] = Aprime[j][i] - s;
        }
        else if (j >= i) {
            double s;

            s = 0;
            for (k = 0; k < i; k++) {
                s += L[j][k] * U[k][i];
            }
            L[j][i] = (Aprime[j][i] - s) / U[i][i];
        }
    }
}

---

else if (j >= i)是不必要的，可以只用 else 代替.但是，我们可以拆分 j循环成两个循环，这样都不需要 if/else :

// fix2.c -- split up j's loop to eliminate if/else inside

double U[n][n];
double L[n][n];
double Aprime[n][n];

for (i = 0; i < n; i++) {
    for (j = 0; j <= i; j++) {
        double s = 0;
        for (k = 0; k < j; k++)
            s += L[j][k] * U[k][i];
        U[j][i] = Aprime[j][i] - s;
    }

    for (; j < n; j++) {
        double s = 0;
        for (k = 0; k < i; k++)
            s += L[j][k] * U[k][i];
        L[j][i] = (Aprime[j][i] - s) / U[i][i];
    }
}

---

U[i][i]在秒中不变 j循环，所以我们可以预先保存它:

// fix3.c -- save off value of U[i][i]

double U[n][n];
double L[n][n];
double Aprime[n][n];

for (i = 0; i < n; i++) {
    for (j = 0; j <= i; j++) {
        double s = 0;
        for (k = 0; k < j; k++)
            s += L[j][k] * U[k][i];
        U[j][i] = Aprime[j][i] - s;
    }

    double Uii = U[i][i];

    for (; j < n; j++) {
        double s = 0;
        for (k = 0; k < i; k++)
            s += L[j][k] * U[k][i];
        L[j][i] = (Aprime[j][i] - s) / Uii;
    }
}

---

对矩阵的访问可能是缓存性能最差的方式。因此，如果可以翻转维度的分配，则可以大大节省内存访问:

// fix4.c -- transpose matrix coordinates to get _much_ better memory/cache
// performance

double U[n][n];
double L[n][n];
double Aprime[n][n];

for (i = 0; i < n; i++) {
    for (j = 0; j <= i; j++) {
        double s = 0;
        for (k = 0; k < j; k++)
            s += L[k][j] * U[i][k];
        U[i][j] = Aprime[i][j] - s;
    }

    double Uii = U[i][i];

    for (; j < n; j++) {
        double s = 0;
        for (k = 0; k < i; k++)
            s += L[k][j] * U[i][k];
        L[i][j] = (Aprime[i][j] - s) / Uii;
    }
}

---

更新:

In the Op's first k-loop its k<j and in the 2nd k<i don't you have to fix that?

是的，我已经修好了。 fix1.c 的变化太难看了，所以我删除了它并将更改应用到 fix2-fix4在哪里很容易做到。

---

更新#2:

These variables are all local to the function.

如果你的意思是它们是函数范围的[没有 static ]，这表示矩阵不能太大，因为除非代码增加堆栈大小，否则它们会被限制在堆栈大小限制内(例如 8 MB)

尽管矩阵看起来是 VLA [因为 n是小写的]，我忽略了这一点。您可能想尝试使用固定维度数组的测试用例，因为我相信它们可能会更快。

此外，如果矩阵是函数范围的，并且想要并行化事物，您可能需要做(例如)#pragma omp shared(Aprime) shared(U) shared(L) .

对缓存的最大拖累是计算 s 的循环.在 fix4 , 我能够访问 U缓存友好，但是 L访问很差。

I'd need to post a whole lot more if I did include the external context

我也猜到了，所以我推测性地进行了矩阵维度交换，不知道还有多少其他代码需要更改。

我创建了一个新版本来更改 L 上的尺寸回到原来的方式，但将交换的版本保留在其他版本上。这为所有矩阵提供了最佳缓存性能。也就是说，大多数矩阵访问的内部循环使得每次迭代都沿着缓存行递增。

事实上，试一试。它可能会将事情改进到不需要并行的程度。我怀疑代码无论如何都是内存限制的，所以并行可能没有多大帮助。

// fix5.c -- further transpose to fix poor performance on s calc loops
//
// flip the U dimensions back to original

double U[n][n];
double L[n][n];
double Aprime[n][n];

double *Up;
double *Lp;
double *Ap;

for (i = 0; i < n; i++) {
    Ap = Aprime[i];
    Up = U[i];

    for (j = 0; j <= i; j++) {
        double s = 0;
        Lp = L[j];
        for (k = 0; k < j; k++)
            s += Lp[k] * Up[k];
        Up[j] = Ap[j] - s;
    }

    double Uii = Up[i];

    for (; j < n; j++) {
        double s = 0;
        Lp = L[j];
        for (k = 0; k < i; k++)
            s += Lp[k] * Up[k];
        Lp[i] = (Ap[j] - s) / Uii;
    }
}

即使您真的需要原始尺寸，根据其他代码，您也可以转置进入并转置回去。这将使其他代码保持不变，但是，如果此代码确实是一个瓶颈，则额外的转置操作可能足够小，值得这样做。

---

更新 #3:

我对所有版本都运行了基准测试。以下是 n 的耗时和相对于原始的比率等于 1037:

orig: 1.780916929 1.000x
fix1: 3.730602026 0.477x
fix2: 1.743769884 1.021x
fix3: 1.765769482 1.009x
fix4: 1.762100697 1.011x
fix5: 0.452481270 3.936x

比率越高越好。

无论如何，这是我能做的极限了。那么，祝你好运......