c - 乘以int64_t数组的最快方法？

Question

我想向量化两个内存对齐数组的乘法。
我没有找到在AVX / AVX2中乘以64 * 64位的任何方法，所以我只是循环展开和AVX2加载/存储。有更快的方法吗？

注意:我不想保存每个乘法的上半部分结果。

void multiply_vex(long *Gi_vec, long q, long *Gj_vec){

    int i;
    __m256i data_j, data_i;
    __uint64_t *ptr_J = (__uint64_t*)&data_j;
    __uint64_t *ptr_I = (__uint64_t*)&data_i;


    for (i=0; i<BASE_VEX_STOP; i+=4) {
        data_i = _mm256_load_si256((__m256i*)&Gi_vec[i]);
        data_j = _mm256_load_si256((__m256i*)&Gj_vec[i]);

        ptr_I[0] -= ptr_J[0] * q;
        ptr_I[1] -= ptr_J[1] * q;
        ptr_I[2] -= ptr_J[2] * q;
        ptr_I[3] -= ptr_J[3] * q;

        _mm256_store_si256((__m256i*)&Gi_vec[i], data_i);
    }


    for (; i<BASE_DIMENSION; i++)
        Gi_vec[i] -= Gj_vec[i] * q;
}

更新:
我将Haswell微体系结构与两个ICC / GCC编译器一起使用。因此，AVX和AVX2都很好。
在乘法循环展开之后，我用C本征-=替换了_mm256_sub_epi64，在此可以加快速度。目前是ptr_J[0] *= q; ...
我使用__uint64_t，但是是错误。正确的数据类型是__int64_t。

Answer 1

您似乎假设long在您的代码中为64位，但随后也使用__uint64_t。在32bit中，x32 ABI，在Windows中，long是32bit类型。您的标题中提到long long，但是您的代码将其忽略。我想知道您的代码是否假设long是32位。

通过使用AVX256负载，您可以完全射击自己的脚，然后将指针混叠到__m256i上以执行标量操作。 gcc只是放弃了，给了您所要求的糟糕代码: vector 加载，然后是一堆extract和insert指令。您的编写方式意味着还必须解压缩和两个 vector ，以便也可以按标量执行sub，而不是使用vpsubq。

Modern x86 CPUs have very fast L1 cache that can handle two operations per clock。 (Haswell和更高版本:每个时钟两个负载和一个存储)。从同一缓存行执行多个标量加载比 vector 加载和拆包更好。 (但是，不完善的uop调度将吞吐量降低到了其中的84％，请参见下文)

gcc 5.3 -O3 -march=haswell (Godbolt compiler explorer)很好地自动矢量化了一个简单的标量实现。 当没有AVX2时，gcc仍然会愚蠢地使用128b vector 自动矢量化:在Haswell上，这实际上是理想标量64位代码速度的1/2。 (请参阅下面的性能分析，但每个 vector 用2个元素代替4个)。

#include <stdint.h>    // why not use this like a normal person?
#define BASE_VEX_STOP 1024
#define BASE_DIMENSION 1028

// restrict lets the compiler know the arrays don't overlap,
// so it doesn't have to generate a scalar fallback case
void multiply_simple(uint64_t *restrict Gi_vec, uint64_t q, const uint64_t *restrict Gj_vec){
    for (intptr_t i=0; i<BASE_DIMENSION; i++)   // gcc doesn't manage to optimize away the sign-extension from 32bit to pointer-size in the scalar epilogue to handle the last less-than-a-vector elements
        Gi_vec[i] -= Gj_vec[i] * q;
}

内循环:

.L4:
    vmovdqu ymm1, YMMWORD PTR [r9+rax]        # MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B]
    add     rcx, 1    # ivtmp.30,
    vpsrlq  ymm0, ymm1, 32      # tmp174, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B],
    vpmuludq        ymm2, ymm1, ymm3        # tmp173, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B], vect_cst_.25
    vpmuludq        ymm0, ymm0, ymm3        # tmp176, tmp174, vect_cst_.25
    vpmuludq        ymm1, ymm4, ymm1        # tmp177, tmp185, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B]
    vpaddq  ymm0, ymm0, ymm1    # tmp176, tmp176, tmp177
    vmovdqa ymm1, YMMWORD PTR [r8+rax]        # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B]
    vpsllq  ymm0, ymm0, 32      # tmp176, tmp176,
    vpaddq  ymm0, ymm2, ymm0    # vect__13.24, tmp173, tmp176
    vpsubq  ymm0, ymm1, ymm0    # vect__14.26, MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__13.24
    vmovdqa YMMWORD PTR [r8+rax], ymm0        # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__14.26
    add     rax, 32   # ivtmp.32,
    cmp     rcx, r10  # ivtmp.30, bnd.14
    jb      .L4 #,

如果需要，可以将其转换回内部函数，但是让编译器自动矢量化将变得容易得多。我没有尝试分析它是否最佳。

如果通常不使用-O3进行编译，则可以在循环之前使用#pragma omp simd(和-fopenmp)。

当然，它不是标量结尾，而是可能。更快地对Gj_vec的最后32B进行未对齐的加载，并将其存储到Gi_vec的最后32B中，可能与循环中的最后一个存储重叠。 (如果数组小于32B，则仍然需要标量后备。)

Haswell的改进的 vector 固有版本

根据我对Z Boson答案的评论。基于Agner Fog's vector class library code。

Agner Fog的版本保存了一条指令，但是通过使用phadd + pshufd(我在其中使用psrlq / paddq / pand)在shuffle端口上出现了瓶颈。

由于您的一个操作数是常量，因此请确保将set1(q)传递为b，而不是a，以便可以悬挂“bswap”洗牌。

// replace hadd -> shuffle (4 uops) with shift/add/and (3 uops)
// The constant takes 2 insns to generate outside a loop.
__m256i mul64_haswell (__m256i a, __m256i b) {
    // instruction does not exist. Split into 32-bit multiplies
    __m256i bswap   = _mm256_shuffle_epi32(b,0xB1);           // swap H<->L
    __m256i prodlh  = _mm256_mullo_epi32(a,bswap);            // 32 bit L*H products

    // or use pshufb instead of psrlq to reduce port0 pressure on Haswell
    __m256i prodlh2 = _mm256_srli_epi64(prodlh, 32);          // 0  , a0Hb0L,          0, a1Hb1L
    __m256i prodlh3 = _mm256_add_epi32(prodlh2, prodlh);      // xxx, a0Lb0H+a0Hb0L, xxx, a1Lb1H+a1Hb1L
    __m256i prodlh4 = _mm256_and_si256(prodlh3, _mm256_set1_epi64x(0x00000000FFFFFFFF)); // zero high halves

    __m256i prodll  = _mm256_mul_epu32(a,b);                  // a0Lb0L,a1Lb1L, 64 bit unsigned products
    __m256i prod    = _mm256_add_epi64(prodll,prodlh4);       // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32
    return  prod;
}

See it on Godbolt。

请注意，这不包括最终的减法，仅包括乘法。

此版本在Haswell上的性能应比gcc的自动矢量化版本好一些。 (例如，每4个周期一个 vector 而不是每5个周期一个 vector ，这会导致port0吞吐量出现瓶颈。我没有考虑其他瓶颈来解决整个问题，因为这是对答案的补充。)

一个AVX1版本(每个 vector 两个元素)会很糟糕，并且可能仍然比64位标量更糟糕。除非您已经将数据包含在 vector 中，并且想要将结果包含在 vector 中，否则不要这样做(将其提取为标量和反数可能不值得)。

对GCC自动矢量化代码(非固有版本)的性能分析

背景:请参阅Agner Fog's insn tables and microarch guide以及x86标签Wiki中的其他链接。

在AVX512之前(见下文)，这可能仅比标量64位代码快一点:imul r64, m64在Intel CPU上的每个时钟的吞吐量为1(在AMD Bulldozer系列中为每4时钟的吞吐量)。 load / imul / sub-with-memory-dest是Intel CPU上的4个融合域uops(具有可微熔的寻址模式，gcc无法使用)。流水线宽度是每个时钟4个融合域uops，因此即使是大的展开也无法在每个时钟一个时钟上发布。如果展开得足够多，我们将成为加载/存储吞吐量的瓶颈。在Haswell上，每个时钟可以进行2次加载和一个存储，但是窃取加载端口的store-address uops将为lower the throughput to about 81/96 = 84% of that, according to Intel's manual。

因此，也许Haswell的最佳方法是先加载并乘以标量(2 uops)，然后乘以vmovq / pinsrq / vinserti128，以便可以使用vpsubq进行减法。这是8个uops来加载＆乘以所有4个标量，这是7个shuffle uops来将数据放入__m256i(2(movq)+ 4(pinsrq是2 uops)+ 1 vinserti128)，还有3个uops进行 vector 加载/ vpsubq / vector商店。因此，每4乘以18个融合域uops(发出4.5个周期)，但是有7个shuffle uops(执行7个周期)。因此，nvm与纯标量相比并不好。

自动向量化代码对四个值的每个 vector 使用8个 vector ALU指令。在Haswell上，这些微指令中的5个(乘法和移位)只能在端口0上运行，因此，无论您如何展开此算法，它都将每5个周期最多获得一个 vector (即，每5/4个周期一个乘积)。

可以用pshufb(端口5)代替移位，以移动数据并以零移位。 (其他随机播放不支持清零，而不是从输入中复制字节，并且输入中没有任何我们可以复制的已知零。)
paddq / psubq可以在Haswell的端口1/5或Skylake的p015上运行。

Skylake在p01上运行pmuludq和立即数 vector 移位，因此理论上它可以管理每个max(5/2，8/3，11/4)= 11/4 = 2.75周期一个 vector 的吞吐量。因此，它限制了融合域uop的总吞吐量(包括2个 vector 负载和1个 vector 存储)。因此，展开一些循环会有所帮助。由于调度不完善而造成的资源冲突可能会将其瓶颈限制在每个时钟少于4个融合域的范围。希望该循环开销可以在端口6上运行，该端口只能处理某些标量运算，包括add和compare-and-branch，而将端口0/1/5留给 vector ALU op，因为它们接近饱和(8/3 = 2.666个时钟)。不过，加载/存储端口远没有达到饱和。

因此，从理论上讲 Skylake可以每2.75个周期管理一个 vector (加上循环开销)，或者每0.7个周期管理一个乘数，而Haswell的最佳选择(理论上，每1.2个周期有一个标量，理论上每个1.25个周期一个)与 vector )。每〜1.2个周期的标量之一可能需要手动调整的asm循环，因为编译器不知道如何对存储使用单寄存器寻址模式，而对负载使用双寄存器寻址模式(dst + (src-dst)和递增dst)。

同样，如果您的数据在L1缓存中不是很热，那么用更少的指令完成工作就可以使前端领先于执行单元，并在需要数据之前就开始加载。硬件预取不会跨越页面行，因此对于大型数组，甚至在较小的数组中， vector 循环实际上可能会胜过标量。

AVX-512DQ引入了64bx64b-> 64b vector 乘法

如果您添加-mavx512dq，gcc可以使用它自动矢量化。

.L4:
    vmovdqu64       zmm0, ZMMWORD PTR [r8+rax]    # vect__11.23, MEM[base: vectp_Gj_vec.22_86, index: ivtmp.32_76, offset: 0B]
    add     rcx, 1    # ivtmp.30,
    vpmullq zmm1, zmm0, zmm2  # vect__13.24, vect__11.23, vect_cst_.25
    vmovdqa64       zmm0, ZMMWORD PTR [r9+rax]    # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B]
    vpsubq  zmm0, zmm0, zmm1    # vect__14.26, MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__13.24
    vmovdqa64       ZMMWORD PTR [r9+rax], zmm0    # MEM[base: vectp_Gi_vec.19_81, index: ivtmp.32_76, offset: 0B], vect__14.26
    add     rax, 64   # ivtmp.32,
    cmp     rcx, r10  # ivtmp.30, bnd.14
    jb      .L4 #,

因此，如果每条时钟以一条流水线的方式传送这些指令，AVX512DQ(expected to be part of Skylake multi-socket Xeon (Purley) in ~2017)将提供比2倍大的加速比(来自更宽的 vector )。

更新:Skylake-AVX512(又名SKL-X或SKL-SP)对xmm，ymm或zmm vector 以每1.5个周期运行一个VPMULLQ。这是3微秒，延迟时间为15分。 (如果这不是测量故障in the AIDA results，则对于zmm版本，可能会有1c的额外延迟)。
vpmullq的速度比您可以用32位块构建的任何东西都要快得多，因此即使当前的CPU没有64位元素 vector 乘法硬件，也非常有必要为此提供一条指令。 (大概他们使用FMA单位的尾数乘法器。)