openMP COLLAPSE 如何在内部工作？

Question

我正在尝试 openMP 并行性，使用 2 个线程将 2 个矩阵相乘。 我了解外循环并行性是如何工作的（即没有“collapse(2)”的工作原理）。

现在，使用折叠。

#pragma omp parallel for collapse(2) num_threads(2)
    for( i = 0; i < m; i++)
        for( j = 0; j < n; j++)
        {
            s = 0;
            for( k = 0; k < p; k++)
                s += A[i][k] * B[k][j];
            C[i][j] = s;
        }

根据我的收集，折叠将循环“折叠”成一个大循环，然后在大循环中使用线程。所以，对于前面的代码，我认为它相当于这样的：

#pragma omp parallel for num_threads(2)
for (ij = 0; ij <n*m; ij++)
{ 
    i= ij/n; 
    j= mod(ij,n);
    s = 0;
    for( k = 0; k < p; k++)
        s += A[i][k] * B[k][j];
    C[i][j] = s;
}

我的问题是：

1. 它是这样工作的吗？我还没有找到任何关于它的解释 “折叠”循环。
2. 如果是，使用它有什么好处？不 它完全像并行性一样在 2 个线程之间划分作业，而没有 崩溃？。如果不是，那么它是如何工作的？

PS：现在我想得更多，如果 n 是奇数，比如 3，如果没有崩溃，一个线程将有 2 次迭代，而另一个只有一个。这会导致线程的工作不均匀，并且效率会降低。 如果我们使用我的崩溃等效项（如果崩溃确实是这样工作的）每个线程将有“1.5”次迭代。如果 n 非常大，那并不重要，不是吗？更不用说，每次都这样做i= ij/n; j= mod(ij,n);，它会降低性能，不是吗？

Answer 1

OpenMP 规范只是说（Version 4.5 第 58 页）：

如果 collapse 子句指定的参数值大于 1，则该子句适用的关联循环的迭代将折叠到一个更大的迭代空间中，然后根据 schedule 子句进行划分。这些关联循环中迭代的顺序执行决定了折叠迭代空间中迭代的顺序。

所以，基本上你的逻辑是正确的，除了你的代码相当于schedule(static,1) collapse(2) 的情况，即迭代块大小为 1。在一般情况下，大多数 OpenMP 运行时的默认调度为schedule(static)，这意味着块大小将（大约）等于迭代次数除以线程数。然后编译器可能会使用一些优化来实现它，例如为外部循环运行一个固定值的部分内部循环，然后使用完整的内部循环进行整数次外部迭代，然后再次进行部分内部循环。

例如下面的代码：

#pragma omp parallel for collapse(2)
for (int i = 0; i < 100; i++)
    for (int j = 0; j < 100; j++)
        a[100*i+j] = i+j;

被 GCC 的 OpenMP 引擎转化为：

<bb 3>:
i = 0;
j = 0;
D.1626 = __builtin_GOMP_loop_static_start (0, 10000, 1, 0, &.istart0.3, &.iend0.4);
if (D.1626 != 0)
  goto <bb 8>;
else
  goto <bb 5>;

<bb 8>:
.iter.1 = .istart0.3;
.iend0.5 = .iend0.4;
.tem.6 = .iter.1;
D.1630 = .tem.6 % 100;
j = (int) D.1630;
.tem.6 = .tem.6 / 100;
D.1631 = .tem.6 % 100;
i = (int) D.1631;

<bb 4>:
D.1632 = i * 100;
D.1633 = D.1632 + j;
D.1634 = (long unsigned int) D.1633;
D.1635 = D.1634 * 4;
D.1636 = .omp_data_i->a;
D.1637 = D.1636 + D.1635;
D.1638 = i + j;
*D.1637 = D.1638;
.iter.1 = .iter.1 + 1;
if (.iter.1 < .iend0.5)
  goto <bb 10>;
else
  goto <bb 9>;

<bb 9>:
D.1639 = __builtin_GOMP_loop_static_next (&.istart0.3, &.iend0.4);
if (D.1639 != 0)
  goto <bb 8>;
else
  goto <bb 5>;

<bb 10>:
j = j + 1;
if (j <= 99)
  goto <bb 4>;
else
  goto <bb 11>;

<bb 11>:
j = 0;
i = i + 1;
goto <bb 4>;

<bb 5>:
__builtin_GOMP_loop_end_nowait ();

<bb 6>:

这是程序抽象语法树的类 C 表示，可能有点难以阅读，但它的作用是，它只使用一次模运算来计算 i 和 j 的初始值基于对 GOMP_loop_static_start() 的调用确定的迭代块的开始 (.istart0.3)。然后它简单地增加i 和j，因为人们期望实现一个循环嵌套，即增加j 直到它达到100，然后将j 重置为0 并增加i。同时，它还保留了当前迭代次数从.iter.1的折叠迭代空间中，基本同时迭代单个折叠循环和两个嵌套循环。

关于线程数不除以迭代次数的情况，OpenMP标准说：

当没有指定chunk_size时，迭代空间被分成大小大致相等的chunk，每个线程最多分配一个chunk。在这种情况下，块的大小是未指定的。

GCC 实现让具有最高 ID 的线程少做一次迭代。其他可能的分发策略在第 61 页的注释中进行了概述。该列表绝不是详尽的。

Answer 2

标准本身并未指定确切的行为。但是，该标准要求内循环对于外循环的每次迭代都具有完全相同的迭代。这允许进行以下转换：

#pragma omp parallel
{
    int iter_total = m * n;
    int iter_per_thread = 1 + (iter_total - 1) / omp_num_threads(); // ceil
    int iter_start = iter_per_thread * omp_get_thread_num();
    int iter_end = min(iter_iter_start + iter_per_thread, iter_total);

    int ij = iter_start;
    for (int i = iter_start / n;; i++) {
        for (int j = iter_start % n; j < n; j++) {
            // normal loop body
            ij++;
            if (ij == iter_end) {
                goto end;
            }
        }
    }
    end:
}

通过浏览反汇编，我相信这类似于 GCC 所做的。它确实避免了每次迭代的除法/取模，但每个内部迭代器需要一个寄存器和一个加法。当然，不同的调度策略会有所不同。

折叠循环确实增加了可以分配给线程的循环迭代次数，从而有助于负载平衡，甚至首先暴露足够的并行工作。