乘以int64_t数组的最快方法？
我想矢量化两个内存对齐数组的乘法。 我没有找到任何方法在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 -= intrisic _mm256_sub_epi64替换-= ，在那里得到一些加速。 目前，它是ptr_J[0] *= q; ... ptr_J[0] *= q; ...
我使用__uint64_t但是错误 。 正确的数据类型是__int64_t 。
您似乎在代码中假设long为__uint64_t ，但随后也使用__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 // 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
内循环：
.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编译，则可以在循环（和-fopenmp ）之前使用#pragma omp simd 。
当然，它不是标量结尾，而是概率。 更快地执行Gj_vec的最后32B的未对齐加载，并存储到Gi_vec的最后32B，可能与循环中的最后一个存储重叠。 （如果arrays小于32B，仍然需要标量回退。）
改进了Haswell的向量内在版本
从我对Z Boson的回答的评论。 基于Agner Fog的矢量类库代码 。
Agner Fog的版本通过使用phadd + pshufd来保存指令但是在shuffle端口上存在瓶颈，我使用psrlq / paddq / pand。
由于其中一个操作数是常量，请确保将set1(q)作为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 vmovq ），然后是vmovq / pinsrq / vinserti128这样你就可以用vpsubq进行减法。 这是8 uops加载和乘以所有4个标量，7个shuffle uops将数据输入__m256i（2（movq）+ 4（pinsrq是2 uops）+ 1 vinserti128），还有3个uop来执行向量加载/ vpsubq / vector商店。 所以这是每4次乘法的18个融合域uops（发布4.5个循环），但是7个shuffle uops（7个循环执行）。 所以nvm，与纯标量相比，这是不好的。
自动向量化代码对四个值的每个向量使用8个向量ALU指令。 在Haswell上，其中5个uop（乘法和移位）只能在端口0上运行，因此无论您如何展开此算法，它每5个周期最多可以实现一个向量（即每5/4个周期一个乘法）。
可以用pshufb （端口5）替换移位以移动数据并以零移位。 （其他shuffle不支持归零而不是从输入复制一个字节，并且输入中没有任何已知的零可以复制。）
paddq / psubq可以在Haswell上的1/5端口上运行，或者在Skylake上运行psubq 。
Skylake在pmuludq运行pmuludq和立即计数向量移位，因此理论上可以管理每个最大（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缓存中不热，则使用较少的指令完成工作会使前端超出执行单元，并在需要数据之前开始加载。 硬件预取不会跨越页面线，因此矢量循环可能会在大型arrays的实践中击败标量，甚至可能用于较小的arrays 。
AVX-512DQ引入了64bx64b-> 64b向量乘法
如果添加-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（ 预计将在2017年成为Skylake多socketsXeon（Purley）的一部分 ）将提供远大于2倍的加速（来自更宽的向量），如果这些指令每个时钟流水线一次。
更新：Skylake-AVX512（又名SKL-X或SKL-SP）每1.5个循环运行一次VPMULLQ，用于xmm，ymm或zmm向量。 这是3 uops，15c延迟。 （对于zmm版本，可能额外增加1c的延迟，如果这不是AIDA结果中的测量误差。）
vpmullq比你用32位块构建的任何东西vpmullq快得多，所以即使当前的CPU没有64位元素的矢量乘法硬件，也非常值得拥有这样的指令。 （据推测，他们使用FMA单位中的尾数乘数。）
如果您对SIMD 64bx64b到64b（更低）操作感兴趣，请参阅Agner Fog的Vector Class Library中的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在循环中是常量的，你不需要在每次迭代时重新计算一些东西，但你的编译器可能无论如何都会想出来。