如何使用 SSE 内在函数针对打包 32x32 优化 C 代码 => 64 位乘法，并将这些结果的一半解包为（伽罗瓦域）

Question

一段时间以来，我一直在努力解决我正在开发的应用程序中网络编码的性能问题（请参阅优化 SSE -code, 提高网络性能coding-encoding 和 OpenCL 分发）。现在我已经非常接近达到可接受的性能了。这是最内层循环的当前状态（其中花费了超过 99% 的执行时间）：

while(elementIterations-- >0)
{   
    unsigned int firstMessageField = *(currentMessageGaloisFieldsArray++);
    unsigned int secondMessageField = *(currentMessageGaloisFieldsArray++);
    __m128i valuesToMultiply = _mm_set_epi32(0, secondMessageField, 0, firstMessageField);
    __m128i mulitpliedHalves = _mm_mul_epu32(valuesToMultiply, fragmentCoefficentVector);
}

您对如何进一步优化它有什么建议吗？我知道如果没有更多的背景就很难做到，但感谢任何帮助！

Answer 1

现在我醒了，这是我的答案：

在您的原始代码中，瓶颈几乎肯定是_mm_set_epi32。这个单一的内在函数被编译成你程序集中的混乱：

633415EC  xor         edi,edi  
633415EE  movd        xmm3,edi  
...
633415F6  xor         ebx,ebx  
633415F8  movd        xmm4,edi  
633415FC  movd        xmm5,ebx  
63341600  movd        xmm0,esi  
...
6334160B  punpckldq   xmm5,xmm3  
6334160F  punpckldq   xmm0,xmm4 
...
63341618  punpckldq   xmm0,xmm5 

这是什么？ 9 条指令？！？！？！ 纯粹的开销...

另一个看起来奇怪的地方是编译器没有合并添加和加载：

movdqa      xmm3,xmmword ptr [ecx-10h]
paddq       xmm0,xmm3

应该合并到：

paddq       xmm0,xmmword ptr [ecx-10h]

我不确定编译器是否去了脑死亡，或者如果它确实有正当理由这样做……无论如何，与 _mm_set_epi32 相比，这只是一件小事。

免责声明：我将从这里开始展示的代码违反了严格别名。但通常需要非标准兼容方法来实现最大性能。

---

解决方案 1：无矢量化

该解决方案假设 allZero 实际上全为零。

该循环实际上比看起来更简单。由于没有太多算术，最好不要向量化：

//  Test Data
unsigned __int32 fragmentCoefficentVector = 1000000000;

__declspec(align(16)) int currentMessageGaloisFieldsArray_[8] = {10,11,12,13,14,15,16,17};
int *currentMessageGaloisFieldsArray = currentMessageGaloisFieldsArray_;

__m128i currentUnModdedGaloisFieldFragments_[8];
__m128i *currentUnModdedGaloisFieldFragments = currentUnModdedGaloisFieldFragments_;
memset(currentUnModdedGaloisFieldFragments,0,8 * sizeof(__m128i));


int elementIterations = 4;

//  The Loop
while (elementIterations > 0){
    elementIterations -= 1;

    //  Default 32 x 32 -> 64-bit multiply code
    unsigned __int64 r0 = currentMessageGaloisFieldsArray[0] * (unsigned __int64)fragmentCoefficentVector;
    unsigned __int64 r1 = currentMessageGaloisFieldsArray[1] * (unsigned __int64)fragmentCoefficentVector;

    //  Use this for Visual Studio. VS doesn't know how to optimize 32 x 32 -> 64-bit multiply
//    unsigned __int64 r0 = __emulu(currentMessageGaloisFieldsArray[0], fragmentCoefficentVector);
//    unsigned __int64 r1 = __emulu(currentMessageGaloisFieldsArray[1], fragmentCoefficentVector);

    ((__int64*)currentUnModdedGaloisFieldFragments)[0] += r0 & 0x00000000ffffffff;
    ((__int64*)currentUnModdedGaloisFieldFragments)[1] += r0 >> 32;
    ((__int64*)currentUnModdedGaloisFieldFragments)[2] += r1 & 0x00000000ffffffff;
    ((__int64*)currentUnModdedGaloisFieldFragments)[3] += r1 >> 32;

    currentMessageGaloisFieldsArray     += 2;
    currentUnModdedGaloisFieldFragments += 2;
}

在 x64 上编译为这个：

$LL4@main:
mov ecx, DWORD PTR [rbx]
mov rax, r11
add r9, 32                  ; 00000020H
add rbx, 8
mul rcx
mov ecx, DWORD PTR [rbx-4]
mov r8, rax
mov rax, r11
mul rcx
mov ecx, r8d
shr r8, 32                  ; 00000020H
add QWORD PTR [r9-48], rcx
add QWORD PTR [r9-40], r8
mov ecx, eax
shr rax, 32                 ; 00000020H
add QWORD PTR [r9-24], rax
add QWORD PTR [r9-32], rcx
dec r10
jne SHORT $LL4@main

在 x86 上编译为这个：

$LL4@main:
mov eax, DWORD PTR [esi]
mul DWORD PTR _fragmentCoefficentVector$[esp+224]
mov ebx, eax
mov eax, DWORD PTR [esi+4]
mov DWORD PTR _r0$31463[esp+228], edx
mul DWORD PTR _fragmentCoefficentVector$[esp+224]
add DWORD PTR [ecx-16], ebx
mov ebx, DWORD PTR _r0$31463[esp+228]
adc DWORD PTR [ecx-12], edi
add DWORD PTR [ecx-8], ebx
adc DWORD PTR [ecx-4], edi
add DWORD PTR [ecx], eax
adc DWORD PTR [ecx+4], edi
add DWORD PTR [ecx+8], edx
adc DWORD PTR [ecx+12], edi
add esi, 8
add ecx, 32                 ; 00000020H
dec DWORD PTR tv150[esp+224]
jne SHORT $LL4@main

这两个都可能已经比原始（SSE）代码更快... 在 x64 上，展开它会让它变得更好。

---

解决方案 2：SSE2 Integer Shuffle

此解决方案将循环展开为 2 次迭代：

//  Test Data
__m128i allZero = _mm_setzero_si128();
__m128i fragmentCoefficentVector = _mm_set1_epi32(1000000000);

__declspec(align(16)) int currentMessageGaloisFieldsArray_[8] = {10,11,12,13,14,15,16,17};
int *currentMessageGaloisFieldsArray = currentMessageGaloisFieldsArray_;

__m128i currentUnModdedGaloisFieldFragments_[8];
__m128i *currentUnModdedGaloisFieldFragments = currentUnModdedGaloisFieldFragments_;
memset(currentUnModdedGaloisFieldFragments,0,8 * sizeof(__m128i));


int elementIterations = 4;

//  The Loop
while(elementIterations > 1){   
    elementIterations -= 2;

    //  Load 4 elements. If needed use unaligned load instead.
    //      messageField = {a, b, c, d}
    __m128i messageField = _mm_load_si128((__m128i*)currentMessageGaloisFieldsArray);

    //  Get into this form:
    //      values0 = {a, x, b, x}
    //      values1 = {c, x, d, x}
    __m128i values0 = _mm_shuffle_epi32(messageField,216);
    __m128i values1 = _mm_shuffle_epi32(messageField,114);

    //  Multiply by "fragmentCoefficentVector"
    values0 = _mm_mul_epu32(values0, fragmentCoefficentVector);
    values1 = _mm_mul_epu32(values1, fragmentCoefficentVector);

    __m128i halves0 = _mm_unpacklo_epi32(values0, allZero);
    __m128i halves1 = _mm_unpackhi_epi32(values0, allZero);
    __m128i halves2 = _mm_unpacklo_epi32(values1, allZero);
    __m128i halves3 = _mm_unpackhi_epi32(values1, allZero);


    halves0 = _mm_add_epi64(halves0, currentUnModdedGaloisFieldFragments[0]);
    halves1 = _mm_add_epi64(halves1, currentUnModdedGaloisFieldFragments[1]);
    halves2 = _mm_add_epi64(halves2, currentUnModdedGaloisFieldFragments[2]);
    halves3 = _mm_add_epi64(halves3, currentUnModdedGaloisFieldFragments[3]);

    currentUnModdedGaloisFieldFragments[0] = halves0;
    currentUnModdedGaloisFieldFragments[1] = halves1;
    currentUnModdedGaloisFieldFragments[2] = halves2;
    currentUnModdedGaloisFieldFragments[3] = halves3;

    currentMessageGaloisFieldsArray     += 4;
    currentUnModdedGaloisFieldFragments += 4;
}

编译为此 (x86)：（x64 并没有太大不同）

$LL4@main:
movdqa    xmm1, XMMWORD PTR [esi]
pshufd    xmm0, xmm1, 216               ; 000000d8H
pmuludq   xmm0, xmm3
movdqa    xmm4, xmm0
punpckhdq xmm0, xmm2
paddq     xmm0, XMMWORD PTR [eax-16]
pshufd    xmm1, xmm1, 114               ; 00000072H
movdqa    XMMWORD PTR [eax-16], xmm0
pmuludq   xmm1, xmm3
movdqa    xmm0, xmm1
punpckldq xmm4, xmm2
paddq     xmm4, XMMWORD PTR [eax-32]
punpckldq xmm0, xmm2
paddq     xmm0, XMMWORD PTR [eax]
punpckhdq xmm1, xmm2
paddq     xmm1, XMMWORD PTR [eax+16]
movdqa    XMMWORD PTR [eax-32], xmm4
movdqa    XMMWORD PTR [eax], xmm0
movdqa    XMMWORD PTR [eax+16], xmm1
add       esi, 16                   ; 00000010H
add       eax, 64                   ; 00000040H
dec       ecx
jne       SHORT $LL4@main

仅比两次迭代的非矢量化版本。这使用很少的寄存器，因此即使在 x86 上您也可以进一步展开它。

---

说明：

• 正如 Paul R 提到的，展开两次迭代允许您将初始负载合并到一个 SSE 负载中。这还有一个好处，可以将您的数据存入 SSE 寄存器。
• 由于数据从 SSE 寄存器开始，因此可以用单个 _mm_shuffle_epi32 替换 _mm_set_epi32（在原始代码中编译为大约 9 个指令）。

Answer 2

我建议您将循环展开 2 倍，以便可以使用一个 _mm_load_XXX 加载 4 个 messageField 值，然后将这四个值解压为两个向量对并根据当前循环处理它们。这样，编译器就不会为 _mm_set_epi32 生成大量混乱的代码，并且所有加载和存储都将是 128 位 SSE 加载/存储。这也将使编译器有更多机会在循环内以最佳方式调度指令。