最快的方法来乘以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 inrisic 替换_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指令.你编写它的方式意味着两个向量都必须被解压缩才能sub在标量中进行,而不是使用vpsubq.

现代x86 CPU具有非常快的L1缓存,每个时钟可以处理两个操作.(Haswell及以后:每个时钟有两个负载和一个存储).从同一缓存行执行多个标量加载比向量加载和解包更好.(不完美的uop调度将吞吐量降低到大约84%,但是:见下文)

gcc 5.3 -O3 -march = haswell(Godbolt编译器资源管理器)很好地自动矢量化一个简单的标量实现. 当AVX2不可用时,gcc愚蠢地仍然使用128b向量自动向量化:在Haswell上,这实际上大约是理想标量64位代码速度的1/2. (参见下面的perf分析,但每个向量代替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的矢量类库代码.

Agner Fog的版本通过使用phadd + pshufd来保存指令但是在shuffle端口上存在瓶颈,我使用psrlq/paddq/pand.

由于其中一个操作数是常量,因此请确保传递set1(q)as b,而不是a,因此可以提升"bswap"shuffle.

// 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;
}

在Godbolt看到它.

请注意,这不包括最终减法,只包括乘法.

这个版本在Haswell上的表现应该比gcc的autovectorized版本好一点.(就像每4个周期可能有一个向量而不是每5个周期一个向量,在port0吞吐量方面存在瓶颈.我没有考虑完整问题的其他瓶颈,因为这是答案的后期补充.)

一个AVX1版本(每个向量两个元素)会很糟糕,可能仍然比64位标量差.不要这样做,除非你已经在向量中有你的数据,并希望结果在向量中(提取到标量和返回可能不值得).

GCC自动向量代码的Perf分析(不是内在版本)

背景:参见Agner Fog的insn表和微指南指南,以及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个融合域uop,因此即使是大量的展开也无法在每个时钟发出一次.通过足够的展开,我们将在加载/存储吞吐量方面遇到瓶颈.根据英特尔的手册,Haswell可以实现2个负载和每个时钟一个存储,但是存储地址uops窃取负载端口会将吞吐量降低到大约81/96 = 84%.

所以也许Haswell的最佳方式是加载和乘以标量,(2 uops),然后vmovq/ pinsrq/ vinserti128所以你可以用a减去vpsubq.这是8 uops加载和乘以所有4个标量,7个shuffle uop将数据输入__m256i(2(movq)+ 4(pinsrq是2 uops)+ 1 vinserti128),还有3个uop来执行向量加载/ vpsubq/vector商店.因此,每4次乘法(发布4.5个周期)就有18个融合域uop,但是7个shuffle uops(执行7个周期).所以nvm,与纯标量相比,这并不好.

自动向量化代码对四个值的每个向量使用8个向量ALU指令.在Haswell上,其中5个uop(乘法和移位)只能在端口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个融合域uop.循环开销有望在端口6上运行,端口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(预计将在2017年成为Skylake多插座Xeon(Purley)的一部分)将提供比2x加速(从更宽的矢量)大得多,如果这些指令每个时钟流水线一次.

更新:Skylake-AVX512(又名SKL-X或SKL-SP)每1.5个循环运行一次VPMULLQ,用于xmm,ymm或zmm向量.这是3 uops,15c延迟.(对于zmm版本,可能额外增加1c的延迟,如果这不是AIDA结果中的测量误差.)

vpmullq比使用32位块构建的任何东西要快得多,所以即使当前的CPU没有64位元素的矢量乘法硬件,也非常值得拥有这样的指令.(据推测,他们使用FMA单位中的尾数乘数.)

Answer 2

如果您对 SIMD 64bx64b 到 64b（较低）操作感兴趣，请参阅 Agner Fog 的矢量类库中的 AVX 和 AVX2 解决方案。我将使用数组测试它们，看看它与 GCC 使用通用循环（例如 Peter Cordes 的答案中的循环）所做的比较如何。

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循环内是恒定的，因此您不需要每次迭代都重新计算某些内容，但您的编译器可能会弄清楚这一点。