我想加快我的代码的一部分，但我不认为有一个可能的更好的方法来做以下计算:
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);
}
numBins通常为10，但这位代码被调用的频率很高(每秒80帧的频率，每帧至少调用8次)
我尝试使用一些SSE方法，但它只是稍微加快了这段代码。我认为我可以避免计算两次中点，但我不确定如何。有更好的方法来计算均值和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;
}
我可能做错了什么，因为我没有很多的改进，因为我期待。SSE代码中有什么可能会减慢这个过程吗?
我刚刚意识到您的数据数组开始时是一个int数组，因为您的代码中没有声明。我可以在SSE版本中看到，您从整数开始，然后只存储它的浮点版本。
保持所有的整数将让我们用一个简单的ivec = _mm_add_epi32(ivec, _mm_set1_epi32(4));做循环反向量Aki Suihkonen的答案有一些转换，应该让它优化得更好。特别是，即使没有-ffast-math，自动矢量化器也应该能够做得更多。事实上，它做得很好。你可以用intrinsic做得更好，特别是节省一些32位向量乘法和缩短依赖链。
我的旧答案，基于只是尝试优化您的代码编写，假设FP输入:
您可以使用链接到的算法@Jason，将所有3个循环组合为一个。不过，它可能不会盈利，因为它涉及到一个部门。对于少量的箱子，可能只是循环多次。
首先阅读http://agner.org/optimize/上的指南。他的优化汇编指南中的一些技术将加速您的SSE尝试(我为您编辑了这个问题)。
在可能的情况下组合循环，这样每次加载/存储数据时，您可以对数据做更多的处理。
多个累加器来隐藏循环携带的依赖链的延迟。(在最近的英特尔cpu上，即使是FP添加也需要3个周期。)
- 代替每次迭代的int->float转换，使用float循环计数器和int循环计数器。(每次迭代添加一个
_mm_set1_ps(4.0f)向量。)如果可能的话，在循环中应该避免使用带有可变参数的_mm_set...。它需要几个指令(特别是当每个参数到setr必须单独计算时)
gcc -O3设法自动向量化第一个循环，但不是其他的。使用-O3 -ffast-math，它可以自动矢量化更多。-ffast-math允许它以不同于代码指定的顺序执行FP操作。例如，将vector数组的4个元素相加，只在最后将4个累加器组合。
告诉gcc输入指针按16对齐，可以让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;
}
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 #
因此，gcc不是每次迭代都跳到顶部，而是决定跳到前面复制一个寄存器，然后无条件地将jmp返回到循环的顶部。上循环缓冲区可能会消除这种愚蠢的前端开销，但是gcc应该对循环进行结构化，这样它就不会在每次迭代时复制xmm5->xmm3，然后再复制xmm3->xmm5，因为这太愚蠢了。它应该有条件跳转到循环的顶部。
还请注意gcc用于获取浮点版本循环计数器的技术:从1 2 3 4的整数向量开始，然后添加set1_epi32(4)。将其用作已打包的int->float cvtdq2ps的输入。在英特尔硬件上，该指令运行在FP-add端口上，并且具有3个周期的延迟，与封装FP add. gcc问题相同。如果只添加一个set1_ps(4.0)的向量会做得更好，尽管这会创建一个3循环的循环依赖链，而不是1循环的向量int add，每次迭代都有3循环的转换分叉。
小迭代计数
你说这通常会恰好用于10个箱子?只有10个箱子的专用版本可以通过避免所有循环开销并将所有内容保存在寄存器中来提供很大的加速。
对于这样小的问题大小，您可以将FP权重放在内存中，而不是每次都使用整数->浮点转换重新计算它们。
另外，10个箱子意味着相对于垂直操作的数量有很多水平操作，因为你只有2.5个向量值的数据。
如果正好10是很常见的，就专门为它指定一个版本。如果16岁以下的孩子很常见，那就专门为他们设计一个版本。(它们可以并且应该共享const float weights[] = { 0.0f, 1.0f, 2.0f, ...};数组。)
您可能希望在专门的小问题版本中使用内在函数，而不是自动向量化。
在您的专用版本中，在数组中有用数据的末尾使用零填充可能仍然是一个好主意。但是，您可以加载最后2个浮点数，并使用movq指令清除向量寄存器的上部64b。( __m128i _mm_cvtsi64_si128 (__int64 a))。将其转换为__m128，就可以了。
正如peterchen所提到的，这些操作对于当前的桌面处理器来说是非常微不足道的。函数是线性的，即O(n)。numBins的典型尺寸是多少?如果它相当大(比如超过1000000)，并行化将会有所帮助。使用像OpenMP这样的库，这可能很简单。如果numBins开始接近MAXINT，你可以考虑GPGPU作为一个选项(CUDA/OpenCL)。
考虑到这一切，您应该尝试分析您的应用程序。如果存在性能约束，很有可能不在这个方法中。Michael Abrash对"高性能代码"的定义极大地帮助了我决定是否/何时进行优化:
在我们能够创建高性能代码之前，我们必须了解什么是高性能。创建高性能软件的目标(并不总是达到的)是使软件能够如此迅速地执行其指定的任务，以至于它能够即时响应，就用户而言。换句话说，高性能代码在理想情况下应该运行得如此之快，以至于代码中的任何进一步改进都是毫无意义的。请注意，上面的定义最强调的是没有提到任何关于使软件尽可能快的内容。
参考:图形编程黑皮书
要计算的整体函数为
std = sqrt(SUM_i { hist[i]/sum * (midpoint_i - mean_midpoint)^2 })
使用标识符
Var (aX + b) = Var (X) * a^2
可以大大降低整体操作的复杂性
1) bin的中点不需要偏移量b
2)不需要预先缩放bin数组元素与bin宽度
和
3)不需要用sum
的倒数来规范化直方图条目优化后的计算如下
float calcVariance(int histBin[], float binWidth)
{
int i;
int sum = 0;
int mid = 0;
int var = 0;
for (i = 0; i < 10; i++)
{
sum += histBin[i];
mid += i*histBin[i];
}
float inv_sum = 1.0f / (float)sum;
float mid_sum = mid * inv_sum;
for (i = 0; i < 10; i++)
{
int diff = i * sum - mid; // because mid is prescaled by sum
var += histBin[i] * diff * diff;
}
return sqrt(float(var) / (float)(sum * sum * sum)) * binWidth;
}
如果float histBin[];
我还将histBin的大小填充为4的倍数，以获得更好的矢量化。
编辑
另一种计算方法是在内循环中使用浮点数:
float inv_sum = 1.0f / (float)sum;
float mid_sum = mid * inv_sum;
float var = 0.0f;
for (i = 0; i < 10; i++)
{
float diff = (float)i - mid_sum;
var += (float)histBin[i] * diff * diff;
}
return sqrt(var * inv_sum) * binWidth;
仅对全局结果进行缩放，并尽可能长时间保持整数
使用Σ(X-m)²/N = ΣX²/N - m²将所有计算分组在一个循环中。
// Accumulate the histogram
int mean= 0, var= 0;
for (int i = 0; i < numBins; ++i)
{
mean+= i * histVec[i];
var+= i * i * histVec[i];
}
// Compute the reduced mean and variance
float fmean= (float(mean) / sum);
float fvar= float(var) / sum - fmean * fmean;
// Rescale
fmean= fmean * binSize + binOffset;
fvar= fvar * binSize * binSize;
所需的整数类型将取决于箱子中的最大值。循环的SSE优化可以利用_mm_madd_epi16指令。
如果箱子的数量小到10，则考虑完全展开循环。预先计算表中的i和i²向量
在幸运的情况下，数据适合16位，和适合32位，累积是用类似
的方式完成的 static short I[16]= { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 0, 0, 0, 0, 0 };
static short I2[16]= { 0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 0, 0, 0, 0, 0, 0 };
// First group
__m128i i= _mm_load_si128((__m128i*)&I[0]);
__m128i i2= _mm_load_si128((__m128i*)&I2[0]);
__m128i h= _mm_load_si128((__m128i*)&inHist[0]);
__m128i mean= _mm_madd_epi16(i, h);
__m128i var= _mm_madd_epi16(i2, h);
// Second group
i= _mm_load_si128((__m128i*)&I[8]);
i2= _mm_load_si128((__m128i*)&I2[8]);
h= _mm_load_si128((__m128i*)&inHist[8]);
mean= _mm_add_epi32(mean, _mm_madd_epi16(i, h));
var= _mm_add_epi32(var, _mm_madd_epi16(i2, h));
注意:未经检查的