10加仑的快速加权平均值和方差

Question

10加仑的快速加权平均值和方差

tre*_*ris 5 c++ optimization sse variance

我想加快我的部分代码,但我认为没有更好的方法来进行以下计算:

float invSum = 1.0f / float(sum);

for (int i = 0; i < numBins; ++i)
{
    histVec[i] *= invSum;
}

for (int i = 0; i < numBins; ++i)
{
    float midPoint = (float)i*binSize + binOffset;
    float f = histVec[i];
    fmean += f * midPoint;
}

for (int i = 0; i < numBins; ++i)
{
    float midPoint = (float)i*binSize + binOffset;
    float f = histVec[i];
    float diff = midPoint - fmean;
    var += f * hwk::sqr(diff);
}

Run Code Online (Sandbox Code Playgroud)

numBins 在for循环中通常为10,但这段代码经常被调用(频率为每秒80帧,每帧至少调用8次)

我尝试使用一些SSE方法,但它只是略微加快了这段代码.我想我可以避免计算两次midPoint,但我不知道如何.有没有更好的方法来计算fmean和var？

这是SSE代码:

// make hist contain a multiple of 4 valid values
    for (int i = numBins; i < ((numBins + 3) & ~3); i++) 
        hist[i] = 0;

// find sum of bins in inHist
    __m128i iSum4 = _mm_set1_epi32(0);
    for (int i = 0; i < numBins; i += 4) 
    {
        __m128i a = *((__m128i *) &inHist[i]);
        iSum4 = _mm_add_epi32(iSum4, a);
    }

    int iSum = iSum4.m128i_i32[0] + iSum4.m128i_i32[1] + iSum4.m128i_i32[2] + iSum4.m128i_i32[3];

    //float stdevB, meanB;

    if (iSum == 0.0f)
    {
        stdev = 0.0;
        mean = 0.0;
    }
    else
    {   
        // Set histVec to normalised values in inHist
        __m128 invSum = _mm_set1_ps(1.0f / float(iSum));
        for (int i = 0; i < numBins; i += 4) 
        {
            __m128i a = *((__m128i *) &inHist[i]);
            __m128  b = _mm_cvtepi32_ps(a);
            __m128  c = _mm_mul_ps(b, invSum);
            _mm_store_ps(&histVec[i], c);
        }

        float binSize = 256.0f / (float)numBins;
        float halfBinSize = binSize * 0.5f;
        float binOffset = halfBinSize;

        __m128 binSizeMask = _mm_set1_ps(binSize);
        __m128 binOffsetMask = _mm_set1_ps(binOffset);
        __m128 fmean4 = _mm_set1_ps(0.0f);
        for (int i = 0; i < numBins; i += 4)
        {
            __m128i idx4 = _mm_set_epi32(i + 3, i + 2, i + 1, i);
            __m128  idx_m128 = _mm_cvtepi32_ps(idx4);
            __m128  histVec4 = _mm_load_ps(&histVec[i]);
            __m128  midPoint4 = _mm_add_ps(_mm_mul_ps(idx_m128, binSizeMask), binOffsetMask);
            fmean4 = _mm_add_ps(fmean4, _mm_mul_ps(histVec4, midPoint4));
        }
        fmean4 = _mm_hadd_ps(fmean4, fmean4); // 01 23 01 23
        fmean4 = _mm_hadd_ps(fmean4, fmean4); // 0123 0123 0123 0123

        float fmean = fmean4.m128_f32[0]; 

        //fmean4 = _mm_set1_ps(fmean);
        __m128 var4 = _mm_set1_ps(0.0f);
        for (int i = 0; i < numBins; i+=4)
        {
            __m128i idx4 = _mm_set_epi32(i + 3, i + 2, i + 1, i);
            __m128  idx_m128 = _mm_cvtepi32_ps(idx4);
            __m128  histVec4 = _mm_load_ps(&histVec[i]);
            __m128  midPoint4 = _mm_add_ps(_mm_mul_ps(idx_m128, binSizeMask), binOffsetMask);
            __m128  diff4 = _mm_sub_ps(midPoint4, fmean4);
            var4 = _mm_add_ps(var4, _mm_mul_ps(histVec4, _mm_mul_ps(diff4, diff4)));
        }

        var4 = _mm_hadd_ps(var4, var4); // 01 23 01 23
        var4 = _mm_hadd_ps(var4, var4); // 0123 0123 0123 0123
        float var = var4.m128_f32[0]; 

        stdev = sqrt(var);
        mean = fmean;
    }

Run Code Online (Sandbox Code Playgroud)

我可能做错了,因为我没有像我期待的那样有很多改进.SSE代码中是否存在可能会降低进程速度的内容？

(编者注:此问题的SSE部分最初被问为/sf/ask/2228647221/,已作为副本关闭.)

Answer 1

Pet*_*des 5

我只是意识到你的数据数组是以int数组开头的,因为你的代码中没有声明.我可以在SSE版本中看到你以整数开头,并且稍后只存储它的浮点版本.

保持一切整数将让我们用简单的ivec = _mm_add_epi32(ivec, _mm_set1_epi32(4)); Aki Suihkonen的答案做循环计数器向量有一些转换,应该让它更好地优化.特别是,自动矢量化器应该能够做得更多甚至没有-ffast-math.事实上,它确实很好.你可以用内在函数做得更好,尤其是保存一些矢量32位乘法并缩短依赖链.

我的旧答案,基于尝试优化您编写的代码,假设FP输入:

您可以使用@Jason链接的算法将所有3个循环合并为一个循环.但是,这可能不是有利可图的,因为它涉及一个分工.对于少量的垃圾箱,可能只是循环多次.

首先阅读http://agner.org/optimize/上的指南.他的优化装配指南中的一些技巧将加速您的SSE尝试(我为您编辑了这个问题).

在可能的情况下组合您的循环,以便在每次加载/存储数据时对数据执行更多操作.
多个累加器来隐藏循环携带的依赖链的延迟.(即使FP添加在最近的Intel CPU上也需要3个周期.)这不适用于像你的情况那样的真正短阵列.
而不是每次迭代的int-> float转换,使用float循环计数器以及int循环计数器.(添加_mm_set1_ps(4.0f)每个迭代的向量.) _mm_set...带有变量args是可能的循环中要避免的东西.它需要几条指令(特别是当每个arg setr必须单独计算时.)

gcc -O3设法自动向量化第一个循环,但不是其他循环.随着-O3 -ffast-math,它自动矢量化更多. -ffast-math允许它以与代码指定的顺序不同的顺序执行FP操作.例如,在向量的4个元素中将数组相加,并且仅在末尾组合4个累加器.

告诉输入指针与16对齐的gcc让gcc自动向量化,开销更少(未对齐部分没有标量循环).

// return mean
float fpstats(float histVec[], float sum, float binSize, float binOffset, long numBins, float *variance_p)
{
    numBins += 3;
    numBins &= ~3;  // round up to multiple of 4.  This is just a quick hack to make the code fast and simple.
    histVec = (float*)__builtin_assume_aligned(histVec, 16);

    float invSum = 1.0f / float(sum);
    float var = 0, fmean = 0;

    for (int i = 0; i < numBins; ++i)
    {
        histVec[i] *= invSum;
        float midPoint = (float)i*binSize + binOffset;
        float f = histVec[i];
        fmean += f * midPoint;
    }

    for (int i = 0; i < numBins; ++i)
    {
        float midPoint = (float)i*binSize + binOffset;
        float f = histVec[i];
        float diff = midPoint - fmean;
//        var += f * hwk::sqr(diff);
        var += f * (diff * diff);
    }
    *variance_p = var;
    return fmean;
}

Run Code Online (Sandbox Code Playgroud)

gcc为第二个循环生成一些奇怪的代码.

        # broadcasting fmean after the 1st loop
        subss   %xmm0, %xmm2    # fmean, D.2466
        shufps  $0, %xmm2, %xmm2        # vect_cst_.16
.L5: ## top of 2nd loop 
        movdqa  %xmm3, %xmm5    # vect_vec_iv_.8, vect_vec_iv_.8
        cvtdq2ps        %xmm3, %xmm3    # vect_vec_iv_.8, vect__32.9
        movq    %rcx, %rsi      # D.2465, D.2467
        addq    $1, %rcx        #, D.2465
        mulps   %xmm1, %xmm3    # vect_cst_.11, vect__33.10
        salq    $4, %rsi        #, D.2467
        paddd   %xmm7, %xmm5    # vect_cst_.7, vect_vec_iv_.8
        addps   %xmm2, %xmm3    # vect_cst_.16, vect_diff_39.15
        mulps   %xmm3, %xmm3    # vect_diff_39.15, vect_powmult_53.17
        mulps   (%rdi,%rsi), %xmm3      # MEM[base: histVec_10, index: _107, offset: 0B], vect__41.18
        addps   %xmm3, %xmm4    # vect__41.18, vect_var_42.19
        cmpq    %rcx, %rax      # D.2465, bnd.26
        ja      .L8     #,   ### <--- This is insane.
        haddps  %xmm4, %xmm4    # vect_var_42.19, tmp160
        haddps  %xmm4, %xmm4    # tmp160, vect_var_42.21
.L2:
        movss   %xmm4, (%rdx)   # var, *variance_p_44(D)
        ret
        .p2align 4,,10
        .p2align 3
.L8:
        movdqa  %xmm5, %xmm3    # vect_vec_iv_.8, vect_vec_iv_.8
        jmp     .L5     #

Run Code Online (Sandbox Code Playgroud)

因此,gcc不是每次迭代都跳回到顶部,而是决定向前跳转以复制寄存器,然后无条件地jmp返回到循环的顶部.uop循环缓冲区可以消除这种愚蠢的前端开销,但是gcc应该已经构造了循环,因此它不会在每次迭代时复制xmm5-> xmm3然后xmm3-> xmm5,因为这很愚蠢.它应该有条件跳转到循环的顶部.

另请注意gcc用于获取循环计数器的浮点版本的技术:以整数向量开始1 2 3 4,然后添加set1_epi32(4).将其用作打包int-> float的输入cvtdq2ps.在Intel HW上,该指令在FP-add端口上运行,并且具有3个周期延迟,与打包FP add相同.gcc prob.最好只添加一个向量set1_ps(4.0),即使这会创建一个3循环循环的依赖链,而不是1循环向量int add,在每次迭代时都有3个周期的转换分支.

小迭代次数

你说这通常会用在10个箱子上吗？只需10个箱的专用版本就可以通过避免所有循环开销并将所有内容保存在寄存器中来实现大幅加速.

使用这个小问题大小,您可以将FP权重放在内存中,而不是每次都使用integer-> float转换重新计算它们.

此外,相对于垂直操作量,10个箱将意味着许多水平操作,因为您只有2个半矢量值的数据.

如果确实10个是非常常见的,那就专门为它设计一个版本.如果16岁以下是常见的,请专门针对该版本.(他们可以而且应该共享const float weights[] = { 0.0f, 1.0f, 2.0f, ...};阵列.)

您可能希望将内在函数用于专门的小问题版本,而不是自动向量化.

在数组中有用数据结束后进行零填充可能仍然是您的专用版本中的一个好主意.但是,您可以加载最后2个浮点数并使用movq指令清除向量寄存器的高64b .(__m128i _mm_cvtsi64_si128 (__int64 a)).把这个投入__m128,你很高兴.

Answer 2

Bil*_*ose 3

正如 peterchen 提到的，这些操作对于当前的桌面处理器来说非常微不足道。该函数是线性的，即 O(n)。的典型尺寸是多少numBins？如果它相当大（例如，超过 1000000），并行化会有所帮助。使用 OpenMP 等库可以很简单。如果numBins开始接近MAXINT，您可以考虑 GPGPU 作为一个选项（CUDA/OpenCL）。

考虑到所有因素，您应该尝试分析您的应用程序。如果存在性能限制，则很可能不在此方法中。Michael Abrash 对“高性能代码”的定义极大地帮助了我确定是否/何时优化：

在我们创建高性能代码之前，我们必须了解什么是高性能。创建高性能软件的目标（并不总是能够实现）是使软件能够快速执行其指定的任务，以便对用户而言能够即时响应。换句话说，高性能代码理想情况下应该运行得如此快，以至于代码的任何进一步改进都是毫无意义的。请注意，上面的定义最强调并没有提到任何关于使软件尽可能快的内容。

参考：图形编程黑皮书

归档时间：	10 年，2 月前
查看次数：	563 次
最近记录：	9 年前