logo标签列表

这个for循环能否并行化？

c++cparallel-processingopenmp
4

我被分配了一些使用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已完成时执行？
（LU分解？）您可以单独计算上三角和下三角矩阵，是的。但是n有多大？OpenMP的开销非常大，我怀疑在大约n = 1e3数量级以下您不会看到太多效果。 - deamentiaemundi
U，L和Aprime是全局范围还是函数范围？这会影响到openmp。有多线程原子更新的考虑，但openmp有处理这个问题的工具。所以，你能贴出更多的代码吗？此外，你有一些不变量，所以看起来可以通过一些重构使其更可并行化。通过拆分事物，你可以从j循环中移除/移动if。另外，n的值是多少？在我看来，它是可以并行化的。 - Craig Estey
在决定并行化某些内容之前，请先查明代码是否有意义。如果代码中只包含一个字母变量，性能可能是最小的问题。首先解决可读性问题，然后修复任何漏洞——当深入挖掘这种烂代码时，可能会出现几个漏洞。例如，此代码将不存在缓存性能，应将其视为漏洞。 - Lundin
@deamentiaemundi 是的，LU分解；-) n可能高达4096，因此计算需要一些时间。 - user2018675
@user2018675 如果您能提供一些数字（Unix的time足够了），特别是涉及内存访问优化的“之前”和“之后”，那就太好了。我发现我的一个旧版本不能充分利用简单（仅有约10％的）并行化，而且我怀疑错误的内存访问是罪魁祸首。 - deamentiaemundi
1个回答
7
在尝试并行化之前，先尝试简化。例如，可以完全消除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;
    }
}

更新：

在Op的第一个k循环中，它是k<j，而在第二个 k<i 中，你需要修复它吗？

是的，我已经修复了。它对于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

比率越高越好。

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

在 Op 的第一个 k 循环中，它是 k<j，而在第二个循环中是 k<i，你不需要修复吗？ - Surt
@Surt 哎呀，你说得对。我漏掉了那个。我会修复它的。这也是我决定分阶段“展示我的工作”的原因之一 :-) - Craig Estey
@Surt 我刚刚更新了，你看一下。是的，楼主的代码可能需要一些修复。我只是在为楼主插入omp内容做准备时进行简化，因为楼主说由于竞争条件，代码无法并行化。但是，openmp有处理竞争的东西。不知道。而且，无论如何，很好地捕捉到了... - Craig Estey
@CraigEstey 谢谢您的反馈。这些变量都是函数内部的局部变量。我没有发布更多内容，因为代码有点乱，如果我包含外部上下文，我需要发布更多内容。感谢您逐步进行操作，您确实帮助我在这个问题上开阔了思路；-)由于我不是编写此代码的人，所以我过于谨慎了。现在我会回去继续工作。谢谢！ - user2018675
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接