有没有可能并行化这个for循环

Is it possible to parallelize this for loop?

本文关键字：for 循环 并行化 有可能      更新时间：2023-10-16

我得到了一些代码来使用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已经完成时才执行?

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

例如，可以完全消除if。

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

注意:在下面的更新#3中，我做了基准测试，从更新#2开始，缓存友好版本fix5的性能比原始版本高出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;
    }
}

更新:

在Op的第一个k-loop是k<j，在第二个k<i，你不需要修复吗?

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

更新# 2:

这些变量都是函数的局部变量。

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

虽然矩阵看起来是VLAs[因为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

比值越高越好。

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

相关文章：