是否有可能将此for循环并行化？

我得到了一些使用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= 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已经完成才能执行else if ？

Answer 1

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

例如， if可以完全消除。

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

注意：在下面的更新＃3中，我做了基准测试，来自更新＃2的缓存友好版本fix5优于原始版本3.9x。

我已经分阶段清理了一些东西，所以你可以看到代码转换。

有了这个，应该可以成功添加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; } } 

---

更新：

在Op的第一个k循环中，它的k而在第二个k ，你不必修复它吗？

是的，我已经修好了。 对于fix1.c ，这是一个太难看的变化，所以我将其删除并将更改应用到fix2-fix4 ，这很容易做到。

---

更新＃2：

这些变量都是函数的本地变量。

如果你的意思是它们是函数作用域[ 没有 static ]，这说明矩阵不能太大，因为除非代码增加了堆栈大小，否则它们仅限于堆栈大小限制（例如8 MB）

虽然矩阵似乎是VLA [因为n是小写的]，但我忽略了这一点。 您可能想尝试使用固定维度数组的测试用例，因为我相信它们可能更快。

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

缓存的最大阻力是计算s的循环。 在fix4 ，我能够访问U缓存友好，但L访问很差。

如果我确实包含外部上下文，我需要发布更多内容

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

我已经创建了一个新版本，它将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 

更高的比率更好。

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