c - 使用 SSE4 向量化点积计算

Question

我正在尝试使用 SSE4 点积改进此代码，但我很难找到解决方案。该函数获取参数 qi 和 tj，其中包含每个 80 个单元格的 float 数组，然后计算点积。返回值是一个包含四个点积的 vector 。所以我想做的是并行计算 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;
}

Answer 1

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

但是，您可以尝试一件事:(我说尝试是因为我还没有计时。)

目前您在 res 上有一个数据依赖链。现在大多数机器上的 vector 加法是 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);