融合三角形循环以实现并行化，计算子索引。

Question

在并行化中常用的一种技术是将嵌套的for循环合并，例如：

for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) {

to

for(int x=0; x<n*n; x++) {
    int i = x/n; int j = x%n;

我想知道如何将三角形循环合并成这样

for(int i=0; i<n; i++) {
   for(int j=0; j<i+1; j++) {

这个程序有 n*(n+1)/2 次迭代。我们将合并后的迭代称为 x。使用二次公式，我得出了以下结果：

for(int x=0; x<(n*(n+1)/2); x++) {      
    int i  = (-1 + sqrt(1.0+8.0*x))/2;
    int j = x - i*(i+1)/2;

与使用正方形循环不同，这需要使用sqrt函数和从整数到浮点数的转换以及从浮点数到整数的转换。

我想知道是否有更简单或更有效的方法？例如，不需要sqrt函数或从整数到浮点数的转换或从浮点数到整数的转换。

编辑：我不想要依赖于先前或下一个迭代的解决方案。我只想要像int i = funci(x) and int j = funcj(x,i)这样的解决方案。

以下是一些代码展示它是如何工作的：

#include <stdio.h>
#include <math.h>
int main() {
    int n = 5;
    int cnt = 0;
    for(int i=0; i<n; i++) {
        for(int j=0; j<i+1; j++) {
            printf("%d: %d %d\n", cnt++, i,j);      
        }
    } printf("\n");

    int nmax = n*(n+1)/2;
    for(int x=0; x<nmax; x++) {     
        int i  = (-1 + sqrt(1.0+8.0*x))/2;
        int j = x - i*(i+1)/2;
        printf("%d: %d %d\n", x,i,j);
    }       
}

- Z boson

---

为什么？如果是为了性能，在最内层循环中调用sqrt()看起来是一个非常不划算的交换。 - unwind

@unwind，它可以用于融合并行for循环。无论如何，融合一个平方循环需要除法（i=x/n，j=x％n），这在现代CPU上并不比sqrt指令慢太多。但这就是问题的关键所在。我能在没有sqrt函数的情况下做到这一点吗？ - Z boson

sqrt 不是唯一的昂贵函数，还包括与双精度浮点数之间的转换。 - MSalters

请注意使用方式；前两个答案对 i 和 j 进行了顺序更新，无法并行化。 - MSalters

如果您只是想要进行速度优化，并且n比较小，那么您可以计算出所有的迭代元组，将它们放入查找表中，然后在这个元组数组上进行迭代。当然，这会增加代码的大小（因此n需要很小）。 - Emilien

Answer 1

考虑到您正在尝试将三角形融合以并行化，非显而易见的解决方案是选择将x映射为(i,j)的非平凡映射：

j |\ i ->
  | \             ____
| |  \    =>    |\\   |
V |___\         |_\\__|

毕竟，您没有按任何特定顺序处理它们，因此确切的映射关系并不重要。

因此，像矩形一样计算x-> i，j ，但如果i > j ，则{ i = N-i，j = N-j } （镜像Y轴，然后镜像X轴）。

   ____
 |\\   |      |\           |\
 |_\\__|  ==> |_\  __  =>  | \
                  / |      |  \
                 /__|      |___\

- MSalters

---

我认为这里有一个打字错误：“但如果 i > N/2” 应该改为 “i > j”，对吗？ - Massimiliano

谢谢，这就是我想要的答案。这是一个聪明的解决方案。OpenMP有一种方法可以融合嵌套循环。我不知道它是否可以融合三角形循环，但我主要想知道如果我已经为正方形循环做了几次，如何手动完成它。 - Z boson

我终于克服了这个问题，并在 答案 中使用了您的解决方案（请查看答案末尾）。三角形循环与我的问题不完全相同。我的问题是带有对角线的左下三角形，而在另一个问题中，三角形是右上角的，不包括对角线。映射并非完全如您所说，但很接近。它变成了 if(j<=i) { i = n - i - 2; j = n - j - 1; }。你是怎么想出这么聪明的答案的？我真的很印象深刻！我很惊讶它没有更多的投票。 - Z boson

@Zboson：从三角形到矩形的几何变换是一种方便证明三角形面积公式 A = h*w/2 的常见方法。编程实际上就是数学。 - MSalters

+1 for the illustration. 我最近为一个分段（更大）的三角形想出了类似的解决方案。我的解决方案是使用迭代计数器，告诉每个线程要跳过多少。即 j = (i + iter) % blockSize; - ofer.sheffer

当j为奇数时，您需要小心一些，否则您会重复处理一些条目。 - Michael Anderson

Answer 2

最合理的形式当然是第一种形式。

话虽如此，融合形式最好使用条件语句实现：

int i = 0; int j = 0;
for(int x=0; x<(n*(n+1)/2); x++) {
  // ...
  ++j;
  if (j>i)
  {
    j = 0;
    ++i;
  }
}

- MSalters

---

在并行循环中，那样做不会轻松实现。我不想要依赖于前一次或后一次迭代的解决方案。它们必须是独立的。非常抱歉我在问题中没有明确说明这一点。我已经更新了我的问题。 - Z boson

Answer 3

我在想是否有更简单或更有效的方法来完成这个任务？

是的，你一开始的代码就可以。请记住以下几点：

• 不存在浮点运算比普通整数更快的情况。
• 然而，存在许多情况，其中浮点数比普通整数慢得多。无论是否有FPU。
• 在大多数系统上，浮点变量通常比普通整数更大，因此仅出于这个原因就更慢。
• 代码的第一个版本可能对缓存内存最友好。对于任何手动优化的情况，这完全取决于您使用的CPU。
• 除法通常在大多数系统上都很慢，无论是对普通整数还是浮点数进行除法运算。
• 任何形式的复杂算术都比简单计数慢。

因此，对于世界上任何给定的CPU，你的第二个示例基本上保证比第一个示例要慢得多。此外，它也完全不可读。

- Lundin

---

你是否认为在并行化中合并平方循环通常没有帮助？ - Z boson

@Zboson 我的观点是，复杂的代码永远不可能比简单的代码更快。在你对程序进行基准测试并找到瓶颈之前，试图手动优化毫无意义，而且你必须有一个具体的系统和CPU，并且对这个系统/CPU有非常深入的了解。尽管结果表明，从70年代的古老8位MCU到现代64位怪兽的任何CPU上，你所编写的代码都会很慢。 - Lundin

我理解你的所有观点。我可能不应该使用“更有效”的词语。通常，在循环中计算/负载所需的时间比计算迭代器的时间要长得多，因此融合的效率并不那么重要。但是，融合可以用于负载分配。我主要关注与我提出的融合循环不同的解决方案，正如你所说的那样，“完全无法阅读”。如果解决方案更有效，那就更好了。 - Z boson