乘以 int64_t 数组的最快方法？
Fastest way to multiply an array of int64_t?
我想向量化两个内存对齐数组的乘法。
我没有找到任何方法在 AVX/AVX2 中乘以 64*64 位，所以我只是做了循环展开和 AVX2 loads/stores。有更快的方法吗？
注意：我不想保存每次乘法的高半部分
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; ...
我使用 __uint64_t 但是 错误 。正确的数据类型是 __int64_t.
您似乎假设 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)auto-vectorizes一个简单的标量实现相当不错。 当 AVX2 不可用时，gcc 仍然愚蠢地 auto-vectorizes 使用 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）。
当然，它可能不是标量结尾，而是 prob。更快地对 Gj_vec 的最后 32B 进行未对齐加载，并存储到 Gi_vec 的最后 32B，可能与循环中的最后一个存储重叠。 （如果数组小于 32B，仍然需要标量回退。）
改进了 Haswell 的矢量内在版本
来自我对 Z Boson 的回答的评论。基于 Agner Fog's vector class library code.
Agner Fog 的版本通过使用 phadd + pshufd 节省了一条指令，但在 shuffle 端口上存在瓶颈，我使用 psrlq / paddq / pand。
由于您的操作数之一是常量，请确保将 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;
}
请注意，这不包括最后的减法，只包括乘法。
这个版本在 Haswell 上的表现应该比 gcc 的自动矢量化版本好一点。 （比如可能每 4 个周期一个向量而不是每 5 个周期一个向量，在 port0 吞吐量上存在瓶颈。我没有考虑整个问题的其他瓶颈，因为这是对答案的后期补充。）
AVX1 版本（每个向量两个元素）会很糟糕，并且可能仍然比 64 位标量差。不要这样做，除非你已经有向量中的数据，并且想要向量中的结果（提取到标量并返回可能不值得）。
对 GCC 的自动向量化代码（不是内部版本）进行 Perf 分析
背景：参见 Agner Fog's insn tables and microarch guide, and other links in the x86 标签 wiki。
直到 AVX512（见下文），这可能只比标量 64 位代码快一点点：imul r64, m64 在 Intel CPU 上每时钟一个吞吐量（但在 AMD 上每 4 个时钟一个吞吐量 Bulldozer-family). load/imul/sub-with-memory-dest 在 Intel CPU 上是 4 fused-domain 微指令（具有可以 micro-fuse 的寻址模式，gcc 无法使用）。流水线宽度为每个时钟 4 fused-domain 微指令，因此即使是大的展开也无法在 one-per-clock 上发布。展开足够多后，我们将在 load/store 吞吐量上遇到瓶颈。在 Haswell 上每个时钟 2 个加载和一个存储是可能的，但是 stores-address uops 窃取加载端口将 lower the throughput to about 81/96 = 84% of that, according to Intel's manual.
因此，对于 Haswell 来说，最好的方法可能是用标量加载和乘法，(2 微指令)，然后 vmovq / pinsrq / vinserti128 这样你就可以用 a 做减法vpsubq。这是 8 微指令来加载和乘以所有 4 个标量，7 微指令将数据放入 __m256i（2（movq）+ 4（pinsrq 是 2 微指令）+ 1 vinserti128），还有 3 微指令来进行矢量加载/vpsubq/向量存储。所以这是每 4 次乘法 18 fused-domain 微指令（发出 4.5 个周期），但是 7 个随机微指令（执行 7 个周期）。所以nvm，这和纯标量相比是不行的。
自动向量化代码对每个包含四个值的向量使用 8 个向量 ALU 指令。在 Haswell 上，其中 5 个 uops（乘法和移位）只能在端口 0 上 运行，因此无论您如何展开该算法，它最多只能实现每 5 个周期一个向量（即每 5/4 个周期一次乘法） .)
移位可以替换为 pshufb（端口 5）以移动数据并移位为零。 （其他洗牌不支持清零而不是从输入中复制一个字节，并且输入中没有我们可以复制的任何 known-zeros。）
paddq / psubq 可以 运行 在 Haswell 上的端口 1/5 或 Skylake 上的 p015。
天湖运行s pmuludq 和 immediate-count 向量在 p01 上移动，因此理论上它可以管理每个最大（5/2、8/3、11/4）= 一个向量的吞吐量11/4 = 2.75 个周期。因此它在总 fused-domain uop 吞吐量（包括 2 个向量加载和 1 个向量存储）上存在瓶颈。所以一些循环展开会有所帮助。可能来自不完美调度的资源冲突会将其瓶颈限制在每个时钟略低于 4 fused-domain 微指令。循环开销有望在端口 6 上 运行，它只能处理一些标量操作，包括 add 和 compare-and-branch，将端口 0/1/5 留给向量 ALU 操作，因为它们'接近饱和（8/3 = 2.666 个时钟）。不过，load/store 端口远未饱和。
因此，Skylake 理论上可以每 2.75 个周期管理一个向量（加上循环开销），或者每 ~0.7 个周期管理一个乘法，而 Haswell 的最佳选择（每 ~1.2 个周期一个）理论上使用标量的周期，或者理论上每 1.25 个周期使用矢量的周期）。不过，每 ~1.2 个周期的标量可能需要一个 hand-tuned asm 循环，因为编译器不知道如何使用 one-register 存储寻址模式和 two-register 寻址模式对于负载（dst + (src-dst) 和增量 dst）。
此外，如果您的数据在 L1 缓存中不热，则使用更少的指令完成工作可以让前端领先于执行单元，并在数据被加载之前开始加载需要。硬件预取不会跨页行，因此矢量循环在实践中可能会胜过大型数组的标量，甚至可能是较小的数组.
AVX-512DQ 引入了 64bx64b->64b 向量乘法
gcc 可以 auto-vectorize 使用它，如果你添加 -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）运行s VPMULLQ，对于 xmm、ymm 或 zmm 向量，每 1.5 个周期一个。它是 3 微指令，延迟为 15c。 （zmm 版本可能有额外的 1c 延迟，如果这不是测量故障 in the AIDA results。）
vpmullq 比您可以用 32 位块构建的任何东西都快得多，因此即使当前的 CPU 没有 64-bit-element 也非常值得为此提供一条指令vector-multiply 硬件。 （大概他们在 FMA 单元中使用尾数乘数。）
如果您对 SIMD 64bx64b 到 64b（较低）操作感兴趣，请参阅 Agner Fog 的 Vector Class Library. I would test these with arrays and see how it compares to what GCC does with a generic loop such as the one in Peter Cordes' .
中的 AVX 和 AVX2 解决方案
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 在循环中是常量，你不需要每次迭代都重新计算一些东西，但你的编译器可能会计算出来。
我想向量化两个内存对齐数组的乘法。 我没有找到任何方法在 AVX/AVX2 中乘以 64*64 位，所以我只是做了循环展开和 AVX2 loads/stores。有更快的方法吗？
注意：我不想保存每次乘法的高半部分
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; ...
我使用 __uint64_t 但是 错误 。正确的数据类型是 __int64_t.
您似乎假设 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)auto-vectorizes一个简单的标量实现相当不错。 当 AVX2 不可用时，gcc 仍然愚蠢地 auto-vectorizes 使用 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）。
当然，它可能不是标量结尾，而是 prob。更快地对 Gj_vec 的最后 32B 进行未对齐加载，并存储到 Gi_vec 的最后 32B，可能与循环中的最后一个存储重叠。 （如果数组小于 32B，仍然需要标量回退。）
改进了 Haswell 的矢量内在版本
来自我对 Z Boson 的回答的评论。基于 Agner Fog's vector class library code.
Agner Fog 的版本通过使用 phadd + pshufd 节省了一条指令，但在 shuffle 端口上存在瓶颈，我使用 psrlq / paddq / pand。
由于您的操作数之一是常量，请确保将 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;
}
请注意，这不包括最后的减法，只包括乘法。
这个版本在 Haswell 上的表现应该比 gcc 的自动矢量化版本好一点。 （比如可能每 4 个周期一个向量而不是每 5 个周期一个向量，在 port0 吞吐量上存在瓶颈。我没有考虑整个问题的其他瓶颈，因为这是对答案的后期补充。）
AVX1 版本（每个向量两个元素）会很糟糕，并且可能仍然比 64 位标量差。不要这样做，除非你已经有向量中的数据，并且想要向量中的结果（提取到标量并返回可能不值得）。
对 GCC 的自动向量化代码（不是内部版本）进行 Perf 分析
背景：参见 Agner Fog's insn tables and microarch guide, and other links in the x86 标签 wiki。
直到 AVX512（见下文），这可能只比标量 64 位代码快一点点：imul r64, m64 在 Intel CPU 上每时钟一个吞吐量（但在 AMD 上每 4 个时钟一个吞吐量 Bulldozer-family). load/imul/sub-with-memory-dest 在 Intel CPU 上是 4 fused-domain 微指令（具有可以 micro-fuse 的寻址模式，gcc 无法使用）。流水线宽度为每个时钟 4 fused-domain 微指令，因此即使是大的展开也无法在 one-per-clock 上发布。展开足够多后，我们将在 load/store 吞吐量上遇到瓶颈。在 Haswell 上每个时钟 2 个加载和一个存储是可能的，但是 stores-address uops 窃取加载端口将 lower the throughput to about 81/96 = 84% of that, according to Intel's manual.
因此，对于 Haswell 来说，最好的方法可能是用标量加载和乘法，(2 微指令)，然后 vmovq / pinsrq / vinserti128 这样你就可以用 a 做减法vpsubq。这是 8 微指令来加载和乘以所有 4 个标量，7 微指令将数据放入 __m256i（2（movq）+ 4（pinsrq 是 2 微指令）+ 1 vinserti128），还有 3 微指令来进行矢量加载/vpsubq/向量存储。所以这是每 4 次乘法 18 fused-domain 微指令（发出 4.5 个周期），但是 7 个随机微指令（执行 7 个周期）。所以nvm，这和纯标量相比是不行的。
自动向量化代码对每个包含四个值的向量使用 8 个向量 ALU 指令。在 Haswell 上，其中 5 个 uops（乘法和移位）只能在端口 0 上 运行，因此无论您如何展开该算法，它最多只能实现每 5 个周期一个向量（即每 5/4 个周期一次乘法） .)
移位可以替换为 pshufb（端口 5）以移动数据并移位为零。 （其他洗牌不支持清零而不是从输入中复制一个字节，并且输入中没有我们可以复制的任何 known-zeros。）
paddq / psubq 可以 运行 在 Haswell 上的端口 1/5 或 Skylake 上的 p015。
天湖运行s pmuludq 和 immediate-count 向量在 p01 上移动，因此理论上它可以管理每个最大（5/2、8/3、11/4）= 一个向量的吞吐量11/4 = 2.75 个周期。因此它在总 fused-domain uop 吞吐量（包括 2 个向量加载和 1 个向量存储）上存在瓶颈。所以一些循环展开会有所帮助。可能来自不完美调度的资源冲突会将其瓶颈限制在每个时钟略低于 4 fused-domain 微指令。循环开销有望在端口 6 上 运行，它只能处理一些标量操作，包括 add 和 compare-and-branch，将端口 0/1/5 留给向量 ALU 操作，因为它们'接近饱和（8/3 = 2.666 个时钟）。不过，load/store 端口远未饱和。
因此，Skylake 理论上可以每 2.75 个周期管理一个向量（加上循环开销），或者每 ~0.7 个周期管理一个乘法，而 Haswell 的最佳选择（每 ~1.2 个周期一个）理论上使用标量的周期，或者理论上每 1.25 个周期使用矢量的周期）。不过，每 ~1.2 个周期的标量可能需要一个 hand-tuned asm 循环，因为编译器不知道如何使用 one-register 存储寻址模式和 two-register 寻址模式对于负载（dst + (src-dst) 和增量 dst）。
此外，如果您的数据在 L1 缓存中不热，则使用更少的指令完成工作可以让前端领先于执行单元，并在数据被加载之前开始加载需要。硬件预取不会跨页行，因此矢量循环在实践中可能会胜过大型数组的标量，甚至可能是较小的数组.
AVX-512DQ 引入了 64bx64b->64b 向量乘法
gcc 可以 auto-vectorize 使用它，如果你添加 -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）运行s VPMULLQ，对于 xmm、ymm 或 zmm 向量，每 1.5 个周期一个。它是 3 微指令，延迟为 15c。 （zmm 版本可能有额外的 1c 延迟，如果这不是测量故障 in the AIDA results。）
vpmullq 比您可以用 32 位块构建的任何东西都快得多，因此即使当前的 CPU 没有 64-bit-element 也非常值得为此提供一条指令vector-multiply 硬件。 （大概他们在 FMA 单元中使用尾数乘数。）
如果您对 SIMD 64bx64b 到 64b（较低）操作感兴趣，请参阅 Agner Fog 的 Vector Class Library. I would test these with arrays and see how it compares to what GCC does with a generic loop such as the one in 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 在循环中是常量，你不需要每次迭代都重新计算一些东西，但你的编译器可能会计算出来。