乘以 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。在 32 位中，x32 ABI，在 Windows 上，long 是 32 位类型。你的标题提到了long long，但是你的代码忽略了它。我想知道您的代码是否假设 long 是 32 位的。

使用 AVX256 加载，然后将指针混叠到 __m256i 上以执行标量操作，您完全是在自找麻烦。 gcc 只是放弃并为您提供您要求的糟糕代码：矢量加载，然后是一堆 extract 和 insert 指令。你的写法意味着 both 向量也必须被解包以在标量中执行 sub，而不是使用 vpsubq。

Modern x86 CPUs have very fast L1 cache that can handle two operations per clock。 （Haswell 和更高版本：每个时钟两次加载和一次存储）。从同一缓存行执行多个标量加载比向量加载和解包更好。 （不完善的 uop 调度将吞吐量降低到其中的 84% 左右：见下文）

---

gcc 5.3 -O3 -march=haswell (Godbolt compiler explorer) 可以很好地自动矢量化一个简单的标量实现。 当 AVX2 不可用时，gcc 仍然愚蠢地使用 128b 向量自动向量化：在 Haswell 上，这实际上是理想标量 64 位代码速度的 1/2 左右。（参见下面的性能分析，但每个向量替换 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 的改进向量内在版本

来自我对 Z Boson 的回答。基于Agner Fog's vector class library code。

Agner Fog 的版本在我使用 psrlq / paddq / pand 的地方使用 phadd + pshufd 保存了一条指令，但在 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 个周期一个向量而不是每 5 个周期一个向量，端口 0 吞吐量的瓶颈。我没有考虑整个问题的其他瓶颈，因为这是对答案的后期添加。）

AVX1 版本（每个向量两个元素）会很糟糕，并且可能仍然比 64 位标量差。不要这样做，除非您已经将数据保存在向量中，并且希望将结果保存在向量中（提取到标量并返回可能不值得）。

---

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 个融合域微指令（具有可以微融合的寻址模式，gcc 无法使用）。流水线宽度是每个时钟 4 个融合域 uops，因此即使是大展开也无法让这个问题在每个时钟一个。通过足够的展开，我们将成为加载/存储吞吐量的瓶颈。在 Haswell 上，每个时钟可以进行 2 次加载和 1 次存储，但是窃取加载端口的存储地址 uops 将lower the throughput to about 81/96 = 84% of that, according to Intel's manual。

因此，Haswell 的最佳方法可能是加载和乘以标量（2 微秒），然后是 vmovq / pinsrq / vinserti128，这样您就可以使用 vpsubq 进行减法运算。加载和乘以所有 4 个标量需要 8 个 uops，将数据放入 __m256i（2（movq）+4（pinsrq 是 2 个 uops）+1 个 vinserti128）需要 7 个 uop，另外 3 个 uops 可以进行向量加载/vpsubq/向量店铺。因此，每 4 次乘法有 18 个融合域微指令（发出 4.5 个周期），但有 7 个随机微指令（执行 7 个周期）。所以 nvm，这与纯标量相比是不好的。

---

自动向量化代码对每个包含四个值的向量使用 8 个向量 ALU 指令。在 Haswell 上，其中 5 个微指令（乘法和移位）只能在端口 0 上运行，因此无论您如何展开该算法，它最多只能实现每 5 个周期一个向量（即每 5/4 个周期一个乘法）。

可以用pshufb（端口5）替换移位以移动数据并移入零。 （其他 shuffle 不支持归零，而是从输入中复制一个字节，并且输入中没有我们可以复制的任何已知零。）

paddq / psubq 可以在 Haswell 上的端口 1/5 或 Skylake 上的 p015 上运行。

Skylake 运行 pmuludq 并且立即计数向量在 p01 上移动，因此理论上它可以管理每个最大值（5/2、8/3、11/4）= 11/4 = 2.75 的一个向量的吞吐量循环。因此，它成为了总融合域 uop 吞吐量（包括 2 个向量加载和 1 个向量存储）的瓶颈。所以一些循环展开会有所帮助。可能来自不完美调度的资源冲突会将其瓶颈到每个时钟略少于 4 个融合域微指令。循环开销有望在端口 6 上运行，它只能处理一些标量操作，包括 add 和比较和分支，将端口 0/1/5 留给向量 ALU 操作，因为它们接近饱和（8 /3 = 2.666 个时钟）。不过，加载/存储端口远未饱和。

因此，Skylake 理论上可以每 2.75 个周期管理一个向量（加上循环开销），或者每 ~0.7 个周期一个乘法，而 Haswell 的最佳选择（理论上每 ~1.2 个周期一个标量） ，或理论上每 1.25 个周期一个向量）。不过，每 ~1.2 个周期的标量可能需要一个手动调整的 asm 循环，因为编译器不知道如何使用单寄存器寻址模式进行存储，而使用两寄存器寻址模式进行加载（dst + (src-dst)并增加dst)。

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

---

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

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

.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 倍的加速（来自更宽的向量）。

更新：Skylake-AVX512（又名 SKL-X 或 SKL-SP）针对 xmm、ymm 或 zmm 向量以每 1.5 个周期运行一个 VPMULLQ。它是 3 微指令，延迟为 15c。 （zmm 版本可能会有额外 1c 的延迟，如果这不是测量故障in the AIDA results。）

vpmullq 比你可以用 32 位块构建的任何东西都快得多，因此即使当前的 CPU 没有 64 位元素向量乘法硬件，也非常值得为此提供指令。 （大概他们使用 FMA 单元中的尾数乘数。）

Answer 2

如果您对 SIMD 64bx64b 到 64b（更低）的操作感兴趣，请点击 Agner Fog 的Vector Class Library 的 AVX 和 AVX2 解决方案。我会用数组测试这些，看看它与 GCC 对通用循环所做的比较，例如 Peter Cordes 的answer 中的循环。

AVX（使用 SSE - 您仍然可以使用 -mavx 进行编译以获得 vex 编码）。

// vector operator * : multiply element by element
static inline Vec2q operator * (Vec2q const & a, Vec2q const & b) {
#if INSTRSET >= 5   // SSE4.1 supported
    // instruction does not exist. Split into 32-bit multiplies
    __m128i bswap   = _mm_shuffle_epi32(b,0xB1);           // b0H,b0L,b1H,b1L (swap H<->L)
    __m128i prodlh  = _mm_mullo_epi32(a,bswap);            // a0Lb0H,a0Hb0L,a1Lb1H,a1Hb1L, 32 bit L*H products
    __m128i zero    = _mm_setzero_si128();                 // 0
    __m128i prodlh2 = _mm_hadd_epi32(prodlh,zero);         // a0Lb0H+a0Hb0L,a1Lb1H+a1Hb1L,0,0
    __m128i prodlh3 = _mm_shuffle_epi32(prodlh2,0x73);     // 0, a0Lb0H+a0Hb0L, 0, a1Lb1H+a1Hb1L
    __m128i prodll  = _mm_mul_epu32(a,b);                  // a0Lb0L,a1Lb1L, 64 bit unsigned products
    __m128i prod    = _mm_add_epi64(prodll,prodlh3);       // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32
    return  prod;
#else               // SSE2
    int64_t aa[2], bb[2];
    a.store(aa);                                           // split into elements
    b.store(bb);
    return Vec2q(aa[0]*bb[0], aa[1]*bb[1]);                // multiply elements separetely
#endif
}

AVX2

// vector operator * : multiply element by element
static inline Vec4q operator * (Vec4q const & a, Vec4q const & 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
    __m256i zero    = _mm256_setzero_si256();                 // 0
    __m256i prodlh2 = _mm256_hadd_epi32(prodlh,zero);         // a0Lb0H+a0Hb0L,a1Lb1H+a1Hb1L,0,0
    __m256i prodlh3 = _mm256_shuffle_epi32(prodlh2,0x73);     // 0, a0Lb0H+a0Hb0L, 0, a1Lb1H+a1Hb1L
    __m256i prodll  = _mm256_mul_epu32(a,b);                  // a0Lb0L,a1Lb1L, 64 bit unsigned products
    __m256i prod    = _mm256_add_epi64(prodll,prodlh3);       // a0Lb0L+(a0Lb0H+a0Hb0L)<<32, a1Lb1L+(a1Lb1H+a1Hb1L)<<32
    return  prod;
}

这些函数适用于有符号和无符号 64 位整数。在您的情况下，由于 q 在循环中是恒定的，您不需要每次迭代都重新计算一些东西，但您的编译器可能无论如何都会计算出来。