SSE2：将二维数组中的有符号整数与双精度数相乘并将结果相加在 C 中

Posted 2023-02-16

【中文标题】SSE2：将二维数组中的有符号整数与双精度数相乘并将结果相加在 C 中

【英文标题】：SSE2: Multiplying signed integers from a 2d array with doubles and summing the results in C

【发布时间】：2014-05-04 13:02:42

【问题描述】：

我目前正在尝试对以下代码进行矢量化：

velocity[0] = 0.0;
velocity[1] = 0.0;
velocity[2] = 0.0;
for (int i = 0; i < PARAMQ; i++)

    velocity[0] += currentCell[i] * LATTICEVELOCITIES[i][0];
    velocity[1] += currentCell[i] * LATTICEVELOCITIES[i][1];
    velocity[2] += currentCell[i] * LATTICEVELOCITIES[i][2];

其中 LATTICEVELOCITIES 是一个二维整数数组

 static const int32_t LATTICEVELOCITIES[PARAMQ][3] = 0, -1, -1,
                                      -1, 0, -1,
                                      0, 0, -1,
                                      1, 0, -1,
                                      0, 1, -1,
                                      -1, -1, 0,
                                      0, -1, 0,
                                      1, -1, 0,
                                      -1, 0, 0,
                                      0, 0, 0,
                                      1, 0, 0,
                                      -1, 1, 0,
                                      0, 1, 0,
                                      1, 1, 0,
                                      0, -1, 1,
                                      -1, 0, 1,
                                      0, 0, 1,
                                      1, 0, 1,
                                      0, 1, 1
                                     ;

currentCell 是一个双精度数组。

不幸的是，我能找到的所有示例都只处理相同类型的数组，我不知道如何只将两个整数加载到一个 128 位寄存器中，然后将它们转换为双精度数。

提前感谢您的帮助。

【问题讨论】：

您的数据组织对于矢量化来说不是最理想的——您能改变它吗？

我宁愿不改变它，但如果有必要我可以。你有什么建议？

如果您可以添加 LATTICEVELOCITIES 的完整定义，这将有所帮助 - 我假设它类似于 int32_t LATTICEVELOCITIES[PARAMQ][3]; ？

我编辑了原始帖子以现在包含定义。

好的 - 谢谢 - 我们可以将定义更改为：static const double LATTICEVELOCITIES[PARAMQ][4];（第四个元素只是一个虚拟元素），还是更好：static const double LATTICEVELOCITIES[3][PARAMQ];？

【参考方案1】：

首先，根据上面的 cmets，我将假设转置 LATTICEVELOCITIES 是可以的：

static const int32_t LATTICEVELOCITIES[3][PARAMQ] = 
     0, -1, 0, 1, 0, -1, 0, 1, -1, 0, 1, -1, 0, 1, 0, -1, 0, 1, 0 ,
     -1, 0, 0, 0, 1, -1, -1, -1, 0, 0, 0, 1, 1, 1, -1, 0, 0, 0, 1 ,
     -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1 
;

现在让我们遍历数据，每次迭代处理两个元素，最后有一个标量迭代来处理最后一个（奇数）元素：

__m128d v0, v2, v2;
v0 = v1 = v2 = _mm_setzero_pd();
for (int i = 0; i < PARAMQ - 1; i += 2)

    __m128d vc, vl0, vl1, vl2;
    __m128i vtemp;

    vc = _mm_loadu_pd(&currentCell[i]);
    vtemp = _mm_loadu_si128((__m128i *)&LATTICEVELOCITIES[0][i]);
    vl0 = _mm_cvtepi32_pd(vtemp);
    vtemp = _mm_loadu_si128((__m128i *)&LATTICEVELOCITIES[1][i]);
    vl1 = _mm_cvtepi32_pd(vtemp);
    vtemp = _mm_loadu_si128((__m128i *)&LATTICEVELOCITIES[2][i]);
    vl2 = _mm_cvtepi32_pd(vtemp);
    v0 = _mm_add_pd(v0, _mm_mul_pd(vc, vl0));
    v1 = _mm_add_pd(v1, _mm_mul_pd(vc, vl1));
    v2 = _mm_add_pd(v2, _mm_mul_pd(vc, vl2));

v0 = _mm_hadd_pd(v0, v0);
v1 = _mm_hadd_pd(v1, v1);
v2 = _mm_hadd_pd(v2, v2);
_mm_store_sd(&velocity[0], v0);
_mm_store_sd(&velocity[1], v1);
_mm_store_sd(&velocity[2], v2);
if (i < PARAMQ)

    velocity[0] += currentCell[i] * LATTICEVELOCITIES[0][i];
    velocity[1] += currentCell[i] * LATTICEVELOCITIES[1][i];
    velocity[2] += currentCell[i] * LATTICEVELOCITIES[2][i];

请注意，这是完全未经测试的代码 - 对于需要修复的拼写错误或错误，我们深表歉意，但基本思想应该是合理的。

还请注意，您应该针对等效的标量代码对此进行测试和基准测试 - 现代 CPU 通常有两个 FPU，因此从 SSE 中可能没有太多收获。如果您可以假设 Sandy Bridge/Ivy Bridge/Haswell 或更高版本，那么 AVX/AVX2 实现应该会做得更好。

【讨论】：

非常感谢！它运行良好，节省了 20 秒的计算时间。我还不熟悉 AVX，但我可能会研究一下。