使用 SSE4 向量化点积计算答案

【问题标题】：Vectorizing Dot Product Calculation using SSE4使用 SSE4 向量化点积计算
【发布时间】：2012-07-04 11:49:13
【问题描述】：

我正在尝试使用 SSE4 点积改进此代码，但我很难找到解决方案。此函数获取参数 qi 和 tj，其中包含每个 80 个单元格的浮点数组，然后计算点积。返回值是一个具有四个点积的向量。所以我要做的是并行计算二十个值的四个点积。

你知道如何改进这段代码吗？

inline __m128 ScalarProd20Vec(__m128* qi, __m128* tj)
{
    __m128 res=_mm_add_ps(_mm_mul_ps(tj[0],qi[0]),_mm_mul_ps(tj[1],qi[1]));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[2],qi[2]),_mm_mul_ps(tj[3],qi[3])));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[4],qi[4]),_mm_mul_ps(tj[5],qi[5])));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[6],qi[6]),_mm_mul_ps(tj[7],qi[7])));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[8],qi[8]),_mm_mul_ps(tj[9],qi[9])));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[10],qi[10]),_mm_mul_ps(tj[11],qi[11])));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[12],qi[12]),_mm_mul_ps(tj[13],qi[13])));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[14],qi[14]),_mm_mul_ps(tj[15],qi[15])));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[16],qi[16]),_mm_mul_ps(tj[17],qi[17])));
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[18],qi[18]),_mm_mul_ps(tj[19],qi[19])));
    return res;
}

【问题讨论】：

标签： c performance sse dot-product

【解决方案1】：

在我在 SO 上看到的数百个 SSE 示例中，您的代码是少数从一开始就处于良好状态的代码之一。您不需要 SSE4 点积指令。（你可以做得更好！）

但是，您可以尝试一件事：（我说尝试是因为我还没有计时。）

目前您在res 上有一个数据依赖链。当今大多数机器上的向量加法是 3-4 个周期。因此，您的代码至少需要运行 30 个周期才能运行：

(10 additions on critical path) * (3 cycles addps latency) = 30 cycles

您可以按如下方式对res 变量进行节点拆分：

__m128 res0 = _mm_add_ps(_mm_mul_ps(tj[ 0],qi[ 0]),_mm_mul_ps(tj[ 1],qi[ 1]));
__m128 res1 = _mm_add_ps(_mm_mul_ps(tj[ 2],qi[ 2]),_mm_mul_ps(tj[ 3],qi[ 3]));

res0 = _mm_add_ps(res0,_mm_add_ps(_mm_mul_ps(tj[ 4],qi[ 4]),_mm_mul_ps(tj[ 5],qi[ 5]))); 
res1 = _mm_add_ps(res1,_mm_add_ps(_mm_mul_ps(tj[ 6],qi[ 6]),_mm_mul_ps(tj[ 7],qi[ 7])));

res0 = _mm_add_ps(res0,_mm_add_ps(_mm_mul_ps(tj[ 8],qi[ 8]),_mm_mul_ps(tj[ 9],qi[ 9])));
res1 = _mm_add_ps(res1,_mm_add_ps(_mm_mul_ps(tj[10],qi[10]),_mm_mul_ps(tj[11],qi[11])));

res0 = _mm_add_ps(res0,_mm_add_ps(_mm_mul_ps(tj[12],qi[12]),_mm_mul_ps(tj[13],qi[13])));
res1 = _mm_add_ps(res1,_mm_add_ps(_mm_mul_ps(tj[14],qi[14]),_mm_mul_ps(tj[15],qi[15])));

res0 = _mm_add_ps(res0,_mm_add_ps(_mm_mul_ps(tj[16],qi[16]),_mm_mul_ps(tj[17],qi[17])));
res1 = _mm_add_ps(res1,_mm_add_ps(_mm_mul_ps(tj[18],qi[18]),_mm_mul_ps(tj[19],qi[19])));

return _mm_add_ps(res0,res1);

这几乎将您的关键路径减半。请注意，由于浮点非关联性，这种优化对于编译器来说是非法的。

这是使用 4 路节点拆分和 AMD FMA4 指令的替代版本。如果您不能使用 fused-multiply 添加，请随意拆分它们。它可能仍然比上面的第一个版本更好。

__m128 res0 = _mm_mul_ps(tj[ 0],qi[ 0]);
__m128 res1 = _mm_mul_ps(tj[ 1],qi[ 1]);
__m128 res2 = _mm_mul_ps(tj[ 2],qi[ 2]);
__m128 res3 = _mm_mul_ps(tj[ 3],qi[ 3]);

res0 = _mm_macc_ps(tj[ 4],qi[ 4],res0);
res1 = _mm_macc_ps(tj[ 5],qi[ 5],res1);
res2 = _mm_macc_ps(tj[ 6],qi[ 6],res2);
res3 = _mm_macc_ps(tj[ 7],qi[ 7],res3);

res0 = _mm_macc_ps(tj[ 8],qi[ 8],res0);
res1 = _mm_macc_ps(tj[ 9],qi[ 9],res1);
res2 = _mm_macc_ps(tj[10],qi[10],res2);
res3 = _mm_macc_ps(tj[11],qi[11],res3);

res0 = _mm_macc_ps(tj[12],qi[12],res0);
res1 = _mm_macc_ps(tj[13],qi[13],res1);
res2 = _mm_macc_ps(tj[14],qi[14],res2);
res3 = _mm_macc_ps(tj[15],qi[15],res3);

res0 = _mm_macc_ps(tj[16],qi[16],res0);
res1 = _mm_macc_ps(tj[17],qi[17],res1);
res2 = _mm_macc_ps(tj[18],qi[18],res2);
res3 = _mm_macc_ps(tj[19],qi[19],res3);

res0 = _mm_add_ps(res0,res1);
res2 = _mm_add_ps(res2,res3);

return _mm_add_ps(res0,res2);

【讨论】：

想想吧。有 40 个内存负载。除非您使用的是 Sandy Bridge 处理器，否则您会遇到 40 个周期的瓶颈。所以OP的代码可能已经是最优的了。
关于浮点关联性：编译器标志-ffast-math 经常被低估、误解的害群之马有时会产生奇迹。几乎从人类诞生之日起，AMD 就可以在每个周期执行两次 L1 内存加载，但不幸的是，其他地方的速度都很慢。
非常感谢您的帮助。我的测试结果表明我的代码和你的想法一样快（正如你在评论中提到的那样）。 AMD FMA4 看起来很有趣，但该指令在我的机器上不可用，并且代码必须与 SSE2 兼容。我会用 -ffast-math 试试。
@martins，如果你只能使用 SSE2，为什么要问 SSE4？
好吧，我的说法不够准确。我们有 SSE3 和 SSE4 的标志，但 SSE2 是所需的最小值。

【解决方案2】：

您的数据在内存中的排列格式不适合专门的 SSE4 点积指令 (dpps)。这些指令期望单个向量的维度是相邻的，如下所示：

| dim0 | dim1 | dim2 | ... | dim19 |

而您的数据似乎具有相互交错的向量：

| v0-dim0 | v1-dim0 | v2-dim0 | v3-dim0 | v0-dim1 | ...

您当前的一般方法似乎是合适的 - 您可以通过重新排序指令来改进事情，这样乘法的结果不会在它们生成后立即使用，但实际上编译器应该能够解决这个问题自己的。

【讨论】：

dpps 仅（可能）对长度为 4 的点积有用，而不是在循环内添加跨多个向量的结果。 inside 循环中不需要任何水平的东西。 dpps 解码为多个微指令以随机播放和水平求和，如果您要在循环中使用它，您仍然需要对其标量结果求和。就像在 Skylake 上一样，3p01 + p5。或者在 Ryzen 上，它是 8 微指令，每 4 个时钟吞吐量 1 个。所以它比使用向量累加器要慢得多。

【解决方案3】：

首先，您可以做的最重要的优化是确保您的编译器已打开所有优化设置。

编译器非常聪明，所以如果把它写成一个循环，它很可能会展开它：

__128 res = _mm_setzero();
for (int i = 0; i < 10; i++) {
  res = _mm_add_ps(res, _mm_add_ps(_mm_mul_ps(tj[2*i], qi[2*i]), _mm_mul_ps(tj[2*i+1], qi[2*i+1])));
}
return res;

（使用 GCC，您需要传递 -funroll-loops，然后它会展开它以一次执行 5 次迭代。）

如果循环版本较慢，您也可以定义一个宏并手动展开它，例如：

__128 res = _mm_setzero();

#define STEP(i) res = _mm_add_ps(res, _mm_add_ps(_mm_mul_ps(tj[2*i], qi[2*i]), _mm_mul_ps(tj[2*i+1], qi[2*i+1])))

STEP(0); STEP(1); STEP(2); STEP(3); STEP(4);
STEP(5); STEP(6); STEP(7); STEP(8); STEP(9);

#undef STEP

return res;

您甚至可以运行从 0 到 20 的循环（或对宏版本执行相同操作），即：

__128 res = _mm_setzero();
for (int i = 0; i < 20; i++) {
  res = _mm_add_ps(res, _mm_mul_ps(tj[i], qi[i]));
}
return res;

（对于 GCC 和 -funroll-loops，这将展开为一次执行 10 次迭代，即与上面的一次两次循环相同。）

【讨论】：