使用SSE4向量化点积计算

Question

使用SSE4向量化点积计算

10

我试图使用SSE4点积来改进这段代码，但是我很难找到一个解决方案。这个函数获取参数qi和tj，它们分别包含80个单元的浮点数组，然后计算点积。返回值是一个向量，其中包含四个点积。所以我想要做的就是并行计算20个值的四个点积。

您有任何想法如何改进这段代码吗？

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;
}

- martin s

3个回答

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个迭代，即与上述两个一起操作的循环相同。

- huon

2

您的数据在内存中的排列格式不适合使用专业的SSE4点积指令(dpps)。这些指令期望单个向量的维度是相邻的，就像这样：

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

您的数据似乎是将向量交错在一起：

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

你目前的一般方法看起来是适当的 - 通过重新安排指令顺序，使得乘法生成的结果不会立即被使用，可能会改善情况，但实际上编译器应该能够自己解决这个问题。

- caf

dpps仅在长度为4的点积中有用，不适用于在循环中添加多个向量结果。循环内部不需要任何水平内容。dpps解码为多个uop以进行洗牌和水平求和，如果您要在循环中使用它，则仍然必须对其标量结果进行求和。就像在Skylake上，3p01 + p5。或者在Ryzen上，它是8个uops，每4个时钟周期1个吞吐量。因此，与使用向量累加器相比，速度要慢得多。 - Peter Cordes

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Mysticial · Accepted Answer

在我看过的所有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指令。如果您不能使用融合乘加操作，请随意拆分它们。它可能仍然比上面的第一个版本更好。

__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);