如何使用 SSE 指令集绝对 2 个双精度或 4 个浮点数？ （最高 SSE4）

Posted 2023-02-16

【中文标题】如何使用 SSE 指令集绝对 2 个双精度或 4 个浮点数？ （最高 SSE4）

【英文标题】：How to absolute 2 double or 4 floats using SSE instruction set? (Up to SSE4)

【发布时间】：2011-04-01 02:44:39

【问题描述】：

这是我尝试使用 SSE 加速的示例 C 代码，这两个数组的长度为 3072 个元素，带有双精度数，如果我不需要双精度数，可以将其降为浮点数。

double sum = 0.0;

for(k = 0; k < 3072; k++) 
    sum += fabs(sima[k] - simb[k]);


double fp = (1.0 - (sum / (255.0 * 1024.0 * 3.0)));

无论如何，我目前的问题是如何在 SSE 寄存器中为双精度或浮点数执行 fabs 步骤，以便我可以将整个计算保留在 SSE 寄存器中，使其保持快速，并且我可以通过部分展开来并行化所有步骤这个循环。

这是我找到的一些资源 fabs() asm 或者可能是这个 flipping the sign - SO 但是第二个资源的弱点需要条件检查。

【参考方案1】：

我建议使用按位并带有掩码。正负值表示相同，只有最高位不同，正值为0，负值为1，见double precision number format。您可以使用以下之一：

inline __m128 abs_ps(__m128 x) 
    static const __m128 sign_mask = _mm_set1_ps(-0.f); // -0.f = 1 << 31
    return _mm_andnot_ps(sign_mask, x);


inline __m128d abs_pd(__m128d x) 
    static const __m128d sign_mask = _mm_set1_pd(-0.); // -0. = 1 << 63
    return _mm_andnot_pd(sign_mask, x); // !sign_mask & x

此外，展开循环以中断循环携带的依赖链可能是个好主意。由于这是非负值的总和，因此求和的顺序并不重要：

double norm(const double* sima, const double* simb) 
__m128d* sima_pd = (__m128d*) sima;
__m128d* simb_pd = (__m128d*) simb;

__m128d sum1 = _mm_setzero_pd();
__m128d sum2 = _mm_setzero_pd();
for(int k = 0; k < 3072/2; k+=2) 
    sum1 += abs_pd(_mm_sub_pd(sima_pd[k], simb_pd[k]));
    sum2 += abs_pd(_mm_sub_pd(sima_pd[k+1], simb_pd[k+1]));


__m128d sum = _mm_add_pd(sum1, sum2);
__m128d hsum = _mm_hadd_pd(sum, sum);
return *(double*)&hsum;

通过展开和打破依赖关系（sum1 和 sum2 现在是独立的），您可以让处理器按顺序执行加法运算。由于指令是在现代 CPU 上流水线化的，因此 CPU 可以在前一个指令完成之前开始处理新指令。此外，按位运算在单独的执行单元上执行，CPU 实际上可以在与加法/减法相同的周期内执行它。我建议Agner Fog's optimization manuals。

最后，我不推荐使用 openMP。循环太小，在多个线程之间分配作业的开销可能大于任何潜在的好处。

【讨论】：

当我有时间的时候，我会看看这个，但我想评论一下 OpenMP 方面......我在这个问题中使用/呈现的 sn-p 是another 大得多的循环内的嵌套循环，我在外循环上使用 OpenMP，而不是内循环 :) 但是，如果 OpenMP 用于这个较小的内循环，你会是正确的。

您好，我能够实现这一点，它确实将速度提高了几个百分点，这很好:) 我真的很喜欢 Agner Fog 的手册，它们真的很棒！

【参考方案2】：

-x 和 x 的最大值应该是 abs(x)。代码如下：

x = _mm_max_ps(_mm_sub_ps(_mm_setzero_ps(), x), x)

【讨论】：

是的，当我发布这个问题时我并没有意识到这个技巧，但是它确实有道理:)

【参考方案3】：

可能最简单的方法如下：

__m128d vsum = _mm_set1_pd(0.0);        // init partial sums
for (k = 0; k < 3072; k += 2)

    __m128d va = _mm_load_pd(&sima[k]); // load 2 doubles from sima, simb
    __m128d vb = _mm_load_pd(&simb[k]);
    __m128d vdiff = _mm_sub_pd(va, vb); // calc diff = sima - simb
    __m128d vnegdiff = mm_sub_pd(_mm_set1_pd(0.0), vdiff); // calc neg diff = 0.0 - diff
    __m128d vabsdiff = _mm_max_pd(vdiff, vnegdiff);        // calc abs diff = max(diff, - diff)
    vsum = _mm_add_pd(vsum, vabsdiff);  // accumulate two partial sums

请注意，这可能不会比现代 x86 CPU 上的标量代码快，后者通常有两个 FPU。但是，如果您可以降低到单精度，那么您很可能会获得 2 倍的吞吐量提升。

另请注意，您需要在循环后将 vsum 中的两个部分和合并为一个标量值，但这很简单，对性能不是很重要。

【讨论】：

我处理部分金额的第二部分花了我一点时间来整理，但总体而言，整个解决方案都有效。谢谢！

@Pharaun：感谢您的反馈 - 您是否对 SSE 代码与原始标量代码进行了基准测试，如果是，您看到了多少性能改进？

@Paul 使用单线程实现我敲掉了 20%，但是通过 openMP 使用多线程它只是一个百分比，所以我仍然需要调整它并考虑转向浮动以获得更高的性能。

@Paul，刚刚将它移到浮点实现上，它在普通 C 上的执行时间缩短了约 15% 到 30%，即使启用了 openMP，它似乎也帮助了我内存带宽 + 四重抽浮子带来了巨大的胜利。

@Pharaun：太好了 - 很高兴它有帮助。如果您可以升级到最新的 Sandy Bridge CPU 并使用 AVX，那么您应该会看到更大的改进。 （与 SSE 的 4 x 浮动相比，AVX 执行 8 x 浮动）。