如何使用SSE2向量化距离计算

Question

A和B是长度为N的向量，其中N可以在20到200之间。 我想要计算这些向量之间距离的平方， 即d^2 = ||A-B||^2。

到目前为止，我有：

float* a = ...;
float* b = ...;
float d2 = 0;

for(int k = 0; k < N; ++k)
{
    float d = a[k] - b[k];
    d2 += d * d;
}

看起来表现良好，但是我已经对我的代码进行了分析，发现这是瓶颈（超过50%的时间都花费在这里）。 我正在使用Visual Studio 2012，在Win 7上，使用这些优化选项：/O2 /Oi /Ot /Oy-。 我的理解是VS2012应该自动向量化该循环（使用SSE2）。 然而，如果我在代码中插入#pragma loop(no_vector)，我并没有感受到明显的减速，所以我猜测该循环没有被向量化。编译器通过以下消息确认了这一点：

  info C5002: loop not vectorized due to reason '1105'

我的问题如下：

1. 是否有可能修复这段代码，以使VS2012可以将其矢量化？
2. 如果不行，那么尝试自己矢量化代码是否有意义？
3. 您能推荐一个网站让我学习SSE2编码吗？
4. 在N的某个值以下，矢量化会适得其反吗？
5. 什么是“reason '1105'”？

- Bull

Answer 1

使用 SSE intrinsics 实现这个很简单：

#include "pmmintrin.h"

__m128 vd2 = _mm_set1_ps(0.0f);
float d2 = 0.0f;
int k;

// process 4 elements per iteration
for (k = 0; k < N - 3; k += 4)
{
    __m128 va = _mm_loadu_ps(&a[k]);
    __m128 vb = _mm_loadu_ps(&b[k]);
    __m128 vd = _mm_sub_ps(va, vb);
    vd = _mm_mul_ps(vd, vd);
    vd2 = _mm_add_ps(vd2, vd);
}

// horizontal sum of 4 partial dot products
vd2 = _mm_hadd_ps(vd2, vd2);
vd2 = _mm_hadd_ps(vd2, vd2);
_mm_store_ss(&d2, vd2);

// clean up any remaining elements
for ( ; k < N; ++k)
{
    float d = a[k] - b[k];
    d2 += d * d;
}

请注意，如果您能保证a和b是16字节对齐的，则可以使用_mm_load_ps而不是_mm_loadu_ps，这可能有助于提高性能，特别是在旧的（尼哈兰之前）CPU上。

此外，请注意，对于像这样的循环，在加载数量相对较少的算术指令时，性能可能会受到内存带宽的限制，并且预期的向量化加速可能无法实现。

- Paul R

---

好的回答！就像我自己做的一样。但是我建议@user2151446确保您使用16字节对齐分配a和b。这样，您可以使用_mm_load_ps而不是_mm_loadu_ps。这将进一步优化代码。 - Øystein Schønning-Johansen

谢谢。这个方法很有效，尽管只有在使用/fp:fast选项时才能与我的原始代码一样快。如果我指定了/fp:fast，那么重要的是在第二个循环中加上#pragma loop(no_vector)，以防止VS尝试对其进行矢量化并且失去第一个循环矢量化带来的一半性能提升。 - Bull

@user2151446：如果你的向量很大并且受到内存带宽限制，那么可能无法从SSE中获得太多好处-我已经在答案中添加了一条注释。 - Paul R

@raxman：其实我对于使用 for (i = 0; i < N - 3; i += 4) 和 for (i = 0; i <= N - 4; i += 4) 这两种方式之间感到犹豫——第一种方式更符合通用的 for 循环写法，而第二种方式则更容易适应不同的向量尺寸（例如，将两个 4 的出现改为 8）。 - Paul R

我通常使用 #define ROUND_UP(x, s) (((x)+((s)-1)) & -(s)) 和 #define ROUND_DOWN(x, s) ((x) & -(s))。我使用 ROUND_UP 来分配数组，这样我就不必担心越界问题，而当我需要担心越界问题时，我使用 ROUND_DOWN。但是你的方法我认为也完全可以。我的意思是，对于一个 8 向量，你只需执行 (i = 0; i <= N - 7; i += 8)。 - user2088790

Answer 2

根据MSDN文档，错误代码1105表示编译器无法将代码简化为规范化指令。对于浮点运算，需要指定/fp:fast选项才能启用任何浮点约简。

- masrtis

---

+1 为/fp:fast选项 - 它确实会使代码向量化。谢谢！然而，即使使用/fp:fast，您的代码实际上需要的时间是我的原始代码（没有/fp:fast）的两倍。new[]/delete[]是一个大问题。通过用堆分配的dist替换为栈分配的缓冲区，可以减少一半的时间，将执行时间恢复到起点。两次循环数据也是一个大问题。事实证明，最好的方法是保留我的原始代码并使用/fp:fast来提高2倍的速度。 - Bull

还有感谢您在MSDN上提供的关于1105错误代码的页面，不确定为什么我找不到它（实际上我已经访问过那个页面！）。如果您删除您的代码，我将接受这个答案。我的原始代码使用/fp:fast是我首选的解决方案。 - Bull

@user2151446 我已经删除了代码，但在我的测试案例中，你的示例代码并没有向量化，因此增加了复杂性。我明天会提出一个相关问题。 - masrtis