使用 SSE 矢量化在 OpenMP 中将内部循环与残差计算并行化

Posted 2023-02-16

【中文标题】使用 SSE 矢量化在 OpenMP 中将内部循环与残差计算并行化

【英文标题】：Parallelizing inner loop with residual calculations in OpenMP with SSE vectorization

【发布时间】：2021-03-30 11:08:25

【问题描述】：

我正在尝试并行化一个程序的内部循环，该程序具有循环范围之外的数据依赖性（最小）。我遇到了一个问题，即残差计算发生在内 j 循环范围之外。如果 j 循环中包含“#pragma omp parallel”部分，即使由于 k 值太低而导致循环根本没有运行，代码也会出错。比如说 (1,2,3)。

for (i = 0; i < 10; i++)
  
    #pragma omp parallel for shared(min) private (j, a, b, storer, arr) //
    for (j = 0; j < k-4; j += 4)
    
      mm_a = _mm_load_ps(&x[j]);
      mm_b = _mm_load_ps(&y[j]);
      mm_a = _mm_add_ps(mm_a, mm_b);
      _mm_store_ps(storer, mm_a);

      #pragma omp critical
      
      if (storer[0] < min)
      
        min = storer[0];
      
      if (storer[1] < min)
      
        min = storer[1];
      
      //etc
      
    
    do
    
        #pragma omp critical
        
        if (x[j]+y[j] < min)
        
          min = x[j]+y[j];
            
         
      
     while (j++ < (k - 1));
    round_min = min
  

【问题讨论】：

您能否更具体地说明您遇到了什么错误？这些是编译时错误还是您得到错误的结果？我的第一个想法是你在min 变量上有一个竞争条件，所以显而易见的解决方案是reduction(min:min) 子句，你必须将它添加到parallel for 指令中。

您应该只使用_mm_min_ps 来获得4 个最小值的向量，并在最后将其减少为一个元素（实际上，由于延迟，您需要多个寄存器），而不是所有分支。这可能有重复。

不过，您需要发布一个实际的minimal reproducible example（没有//etc cmets，也没有使用未声明的变量）。外循环实际上是做什么用的？ x 和 y 改变了吗？

【参考方案1】：

基于j 的循环是一个并行循环，因此您不能在循环之后使用j。尤其如此，因为您明确地将 j 设置为 private，因此仅在线程中本地可见，但在并行区域之外不可见。您可以在并行循环之后使用(k-4+3)/4*4 显式计算剩余j 值的位置。

此外，这里有几个要点：

您可能真的不需要自己矢量化代码：您可以使用omp simd reduction。 OpenMP 可以自动为您完成计算残差计算的所有枯燥工作。此外，代码将是可移植的并且更简单。生成的代码也可能比你的更快。但请注意，某些编译器可能无法对代码进行矢量化（GCC 和 ICC 可以，而 Clang 和 MSVC 通常需要一些帮助）。 关键部分 (omp critical) 非常昂贵。在您的情况下，这只会消除与并行部分相关的任何可能的改进。由于缓存行弹跳，代码可能会变慢。 在此处读取_mm_store_ps 写入的数据效率低效，尽管某些编译器（如 GCC）可能能够理解您的代码逻辑并生成更快的实现（提取车道数据）。李> 水平 SIMD 减少效率低下。使用速度更快且在此处易于使用的垂直模式。

考虑到以上几点，这是一个更正的代码：

for (i = 0; i < 10; i++)

    // Assume min is already initialized correctly here

    #pragma omp parallel for simd reduction(min:min) private(j)
    for (j = 0; j < k; ++j)
    
        const float tmp = x[j] + y[j];
        if(tmp < min)
            min = tmp;
    

    // Use min here

上述代码在 GCC/ICC（均使用 -O3 -fopenmp）、Clang（使用 -O3 -fopenmp -ffastmath）和 MSVC（使用 /O2 /fp:precise -openmp:experimental）上的 x86 架构上正确矢量化。

【讨论】：

您也许可以使用 min = min<tmp ? min : tmp 将编译器手持到更好的 asm 中，以确保他们看到它是无分支的，并且（没有 AVX）让它使用更新累加器寄存器的操作数顺序就地，保存movaps (What is the instruction that gives branchless FP min and max on x86?)。即使没有-ffast-math，也可能会得到很好的结果。 OpenMP 缩减可能会使快速数学变得多余

嗯，即使在没有 openMP godbolt.org/z/d5a6Ma9va 的情况下使用 -ffast-math 时，该三元组实际上也击败了 clang 的自动矢量化。 （这 是 是必要的；min() 实现与 NaN 或有符号零都没有关联。）我简化了外部循环，因此它只是将 1024 个浮点元素减少为一次标量，因此我们可以查看该 asm .

使用 OpenMP，clang 只是并行化，而不是矢量化。 （只是minss，没有minps）godbolt.org/z/s1YzqPq46。即使只使用-O3 -fopenmp，没有-ffast-math，GCC 也会进行矢量化，但这不会让 GCC 使用多个累加器展开以隐藏 MINPS 延迟。 :/ 这对于 GCC 的标准自动矢量化是正常的（没有配置文件引导优化 -fprofile-use），但人们可能希望它的 OpenMP 矢量化器更具侵略性。

是的，确实，Clang 的行为很奇怪（MSVC 也是如此）。这就是我编辑答案以用条件替换三元的原因。我不知道为什么 Clang 自动矢量化器无法识别基于最小值的减少。这似乎是一个错误或内部限制。我很想知道更多关于他为什么会这样的行为。

我猜可能对 OpenMP 的支持有限/不太成熟？并行化外循环并让正常的自动矢量化在内循环上工作可能会很好，这取决于它对局部性的影响。尽管我猜想对于缓存阻塞，您确实希望每个核心都重复接触相同的数据子集，而不是让每个核心循环遍历所有数据。