3

我在搞乱 SSE，试图编写一个函数，将单精度浮点数组的所有值相加。我希望它适用于数组的所有长度，而不仅仅是 4 的倍数，就像网络上几乎所有示例中所假设的那样。我想出了这样的事情：

float sse_sum(const float *x, const size_t n)
{
    const size_t
        steps = n / 4,
        rem = n % 4,
        limit = steps * 4;

    __m128 
        v, // vector of current values of x
        sum = _mm_setzero_ps(0.0f); // sum accumulator

    // perform the main part of the addition
    size_t i;
    for (i = 0; i < limit; i+=4)
    {
        v = _mm_load_ps(&x[i]);
        sum = _mm_add_ps(sum, v);
    }

    // add the last 1 - 3 odd items if necessary, based on the remainder value
    switch(rem)
    {
        case 0: 
            // nothing to do if no elements left
            break;
        case 1: 
            // put 1 remaining value into v, initialize remaining 3 to 0.0
            v = _mm_load_ss(&x[i]);
            sum = _mm_add_ps(sum, v);
            break;
        case 2: 
            // set all 4 to zero
            v = _mm_setzero_ps();
            // load remaining 2 values into lower part of v
            v = _mm_loadl_pi(v, (const __m64 *)(&x[i]));
            sum = _mm_add_ps(sum, v);
            break;
        case 3: 
            // put the last one of the remaining 3 values into v, initialize rest to 0.0
            v = _mm_load_ss(&x[i+2]);
            // copy lower part containing one 0.0 and the value into the higher part
            v = _mm_movelh_ps(v,v);
            // load remaining 2 of the 3 values into lower part, overwriting 
            // old contents                         
            v = _mm_loadl_pi(v, (const __m64*)(&x[i]));     
            sum = _mm_add_ps(sum, v);
            break;
    }

    // add up partial results
    sum = _mm_hadd_ps(sum, sum);
    sum = _mm_hadd_ps(sum, sum);
    __declspec(align(16)) float ret;
    /// and store the final value in a float variable
    _mm_store_ss(&ret, sum);
    return ret; 
}

然后我开始怀疑这不是矫枉过正。我的意思是，我陷入了 SIMD 模式，只需要用 SSE 处理尾部。这很有趣，但是将尾部相加并使用常规浮点运算计算结果不是同样好（而且更简单）吗？我在 SSE 做这件事有什么收获吗？

4

3 回答 3

2

正如所承诺的，我做了一些基准测试。为此，我 _aligned_malloc'd 了一个大小为 100k 的浮点数组，用单个值1.123f填充它，并针对它测试了函数。我写了一个简单的求和函数，它只是将结果累加到一个循环中。接下来，我制作了一个简化的 SSE 求和函数变体，其中水平和尾部添加使用常规浮点数完成：

float sseSimpleSum(const float *x, const size_t n)
{
    /* ... Everything as before, but instead of hadd: */

    // do the horizontal sum on "sum", which is a __m128 by hand
    const float *sumf = (const float*)(&sum);
    float ret = sumf[0] + sumf[1] + sumf[2] + sumf[3];

    // add up the tail
    for (; i < n; ++i)
    {
        ret += x[i];
    }

    return ret; 
}

我的性能没有受到任何影响，有时它甚至看起来有点快，但我发现计时器非常不可靠，所以让我们假设简化的变体等同于复杂的变体。但令人惊讶的是，从 SSE 和幼稚浮点求和函数获得的值存在相当大的差异。我怀疑这是由于四舍五入造成的误差累积，所以我写了一个基于Kahan 算法的函数，它给出了正确的结果，虽然比天真的浮点加法慢得多。为了完整起见，我按照这些思路制作了一个基于 SSE Kahan 的函数：

float SubsetTest::sseKahanSum(const float *x, const size_t n)
{
    /* ... init as before... */

    __m128 
        sum = _mm_setzero_ps(), // sum accumulator
        c   = _mm_setzero_ps(), // correction accumulator
        y, t;

    // perform the main part of the addition
    size_t i;
    for (i = 0; i < limit; i+=4)
    {
        y = _mm_sub_ps(_mm_load_ps(&x[i]), c);
        t = _mm_add_ps(sum, y);
        c = _mm_sub_ps(_mm_sub_ps(t, sum), y);
        sum = t;
    }

    /* ... horizontal and tail sum as before... */
}

以下是从 VC++2010 在 Release 模式下获得的基准测试结果，其中显示了总和的获得值、计算时间以及与正确值相关的错误量：

Kahan: value = 112300, time = 1155, error = 0
Float: value = 112328.78125, time = 323, error = 28.78125
SSE: value = 112304.476563, time = 46, error = 4.4765625
Simple SSE: value = 112304.476563, time = 45,错误 = 4.4765625
Kahan SSE：值 = 112300，时间 = 167，错误 = 0

天真的浮点加法的错误量是巨大的！我怀疑非 Kahan SSE 函数更准确，因为它们相当于成对求和，与直接方法相比，它可以提高准确性。Kahan SSE 是准确的，但只比简单的浮动加法快两倍。

于 2013-03-20T02:20:17.877 回答
2

我会查看 Agner Fog 的矢量类。请参阅VectorClass.pdf的“当数据大小不是向量大小的倍数”部分。他列出了五种不同的方法，并讨论了每种方法的优缺点。 http://www.agner.org/optimize/#vectorclass

一般来说，我这样做的方式是从以下链接中获得的。 http://fastcpp.blogspot.no/2011/04/how-to-unroll-loop-in-c.html

#define ROUND_DOWN(x, s) ((x) & ~((s)-1))
 void add(float* result, const float* a, const float* b, int N) {
 int i = 0;
 for(; i < ROUND_DOWN(N, 4); i+=4) {
    __m128 a4 = _mm_loadu_ps(a + i);
    __m128 b4 = _mm_loadu_ps(b + i);
    __m128 sum = _mm_add_ps(a4, b4);
    _mm_storeu_ps(result + i, sum);
  }
  for(; i < N; i++) {
      result[i] = a[i] + b[i];
  }
}
于 2013-03-20T15:29:01.590 回答
1

在这种情况下，除非您能指出一些真正的性能提升，否则它可能是矫枉过正的。如果您正在使用，gcc那么这个关于使用 gcc 4.7 进行自动矢量化的指南可能是一个不错的选择，尽管它显然是gcc具体的，它不像内在函数那么难看。

于 2013-03-20T00:19:15.397 回答