你的代码似乎假设
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及以后的CPU：每个时钟周期可以进行两次加载和一次存储操作）。从同一缓存行进行多个标量加载要比进行向量加载和解包更好。（不过，由于不完美的微操作调度，吞吐量会降低到大约84%：详见下文）
gcc 5.3 -O3 -march=haswell (Godbolt compiler explorer) 很好地自动向量化了一个简单的标量实现。 当AVX2不可用时，gcc仍然愚蠢地使用128位向量进行自动向量化：在Haswell上，这实际上将是理想标量64位代码速度的约1/2。 （请参见下面的性能分析，但将每个向量替换为2个元素而不是4个）。
#include <stdint.h>
#define BASE_VEX_STOP 1024
#define BASE_DIMENSION 1028
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++)
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，则仍然需要标量回退。）
GCC13和Clang17仍然都使用这种方式进行向量化，使用另外两个`vpmuludq`而不是一个`vpmulld`来完成两个32位的交叉乘积，减少了一些洗牌操作。但它们能够优化掉指针别名的`ptr_I[0] -= ptr_J[0] * q;`中的插入/提取操作，因此编译成与简单标量版本相同的汇编循环。
AVX2的改进向量内在版本
根据我对Z Boson回答的评论。基于Agner Fog的向量类库代码，在2023年进行了各种改进，同时修复了一个错误。（8年来，没有人注意到我将交叉乘积添加到完整结果的底部，而不是与它们的位置值相匹配的上半部分；感谢@Petr实际测试了这一点。）
Agner Fog的版本通过使用phadd + pshufd来节省一条指令，但在洗牌端口上成为瓶颈，而我使用psllq / pand / paddd。在Haswell和后来的Intel上，使用vpmulld将两个交叉乘积打包在一起需要2个uops（与两个单独的vpmuludq指令相同，但输入的洗牌较少），但在Zen 3和后来的AMD上，只需要1个uops。
由于您的操作数之一是常量，请确保将
set1(q) 作为
b 而不是
a 传递，这样 "b_swap" 重排操作就可以被提升。
// replace hadd -> shuffle (4 uops) with shift/and/add (3 uops with less shuffle-port pressure)
// The constant takes 2 insns to generate outside a loop.
__m256i mul64_avx2 (__m256i a, __m256i b)
{
// There is no vpmullq until AVX-512. Split into 32-bit multiplies
// Given a and b composed of high<<32 | low 32-bit halves
// a*b = a_low*(u64)b_low + (u64)(a_high*b_low + a_low*b_high)<<32
// the a_high * b_high product isn't needed for non-widening; its place value is entirely outside the low 64 bits.
__m256i b_swap = _mm256_shuffle_epi32(b, _MM_SHUFFLE(2,3, 0,1)); // swap H<->L
__m256i crossprod = _mm256_mullo_epi32(a, b_swap); // 32-bit L*H and H*L cross-products
__m256i prodlh = _mm256_slli_epi64(crossprod, 32); // bring the low half up to the top of each 64-bit chunk
__m256i prodhl = _mm256_and_si256(crossprod, _mm256_set1_epi64x(0xFFFFFFFF00000000)); // isolate the other, also into the high half were it needs to eventually be
__m256i sumcross = _mm256_add_epi32(prodlh, prodhl); // the sum of the cross products, with the low half of each u64 being 0.
__m256i prodll = _mm256_mul_epu32(a,b); // widening 32x32 => 64-bit low x low products
__m256i prod = _mm256_add_epi32(prodll, sumcross); // add the cross products into the high half of the result
return prod;
}
这已经经过测试并且现在可以正常工作，包括对于具有非零高位和低位的值。
在Godbolt上查看。
请注意，这不包括最后一个问题的减法，只有乘法。
这个版本在Haswell上的表现应该比gcc的自动向量化版本要好一些。（可能是每4个周期一个向量，而不是每5个周期一个向量，瓶颈在于端口0的吞吐量。用一个向量常数替换vpsllq（srli_epi64），或者用4个vpslldq并在加法后进行掩码操作，可以将端口0的瓶颈降低到3个周期。）无论哪种方式，在后来的CPU上，如Skylake和Rocket Lake，
https://uica.uops.info/根据执行端口压力估计，仅对乘法（不包括循环或减法）的吞吐瓶颈为2.67，具有完美调度。我选择了移位版本，这样就不需要额外的向量常数，并且让移位/与操作并行运行，以便在这是依赖链的一部分时获得更好的关键路径延迟，与问题中每个元素在加载和存储之间没有太多其他操作不同，因此乱序执行不需要过多努力在迭代之间重叠。
我们可以完全避免使用向量常量，通过执行
sumcross = ((crossprod >> 32) + crossprod) << 32。这仍然需要3个指令，但AND操作变成了移位操作，因此它是另一个与乘法竞争的微操作（在Intel P-cores上）。它的关键路径延迟更差。但它避免了常量，并且使用了更少的临时寄存器（由于相同的原因，它的ILP更差且延迟更高），因此在作为较大循环体的一部分进行乘法运算时可能会有用，否则可能会用完寄存器。
AVX1或SSE4.1版本（每个向量两个元素）不会比64位标量
imul r64, [mem]更快。除非您已经将数据转换为向量，并且希望结果为向量（与提取到标量然后返回相比可能更好）。或者在具有较慢的
imul r64, r64的CPU上，例如Silvermont系列。在这种情况下，
_mm_add_epi32比
_mm_add_epi64更快。
add_epi64 在任何 AVX2 CPU 上都不会比较慢或者占用更多的代码空间，但是在一些早期的 CPU 上，SSE 版本会比较慢。我在可能的情况下使用 add_epi32，因为它不会更慢，并且可能能够节省一些功耗，或者在一些假设的未来 CPU 上具有更低的延迟。根据 https://uops.info/，自 Haswell 以来的 Intel 和自 Zen 1 以来的 AMD 在 Alder Lake E-cores 上运行 vpaddd 与 vpaddq 完全相同。
GCC自动向量化代码的性能分析（非内部版本）
背景：请参阅Agner Fog的指令表和微架构指南，以及x86标签维基中的其他链接。
直到AVX512（见下文），这可能只比标量64位代码稍微快一点：在Intel CPU上，
imul r64, m64的吞吐量为每个时钟周期一次（但在AMD Bulldozer系列上为每4个时钟周期一次）。在Intel CPU上，load/imul/sub-with-memory-dest是4个融合域微操作（使用可以微融合的寻址模式，但gcc未能使用）。流水线宽度为每个时钟周期4个融合域微操作，因此即使大规模展开也无法使其每个时钟周期发出一次。通过足够的展开，我们将在加载/存储吞吐量上受到瓶颈限制。在Haswell上，每个时钟周期可能有2次加载和1次存储，但存储地址微操作窃取加载端口将使吞吐量降低到约81/96 = 84%（根据Intel手册）。
所以也许Haswell的最佳方式是使用标量进行加载和乘法（2个uops），然后使用vmovq / pinsrq / vinserti128进行减法操作（使用vpsubq）。这样一来，加载和乘法操作需要8个uops来加载和乘法4个标量，使用7个洗牌uops将数据放入__m256i寄存器（2个movq + 4个pinsrq（2个uops）+ 1个vinserti128），然后再使用3个uops进行向量加载/ vpsubq / 向量存储。因此，每4个乘法操作需要18个融合域uops（发射需要4.5个周期），但需要7个洗牌uops（执行需要7个周期）。所以，与纯标量相比，这个方法不好。
自动向量化的代码对于每个包含四个值的向量使用了8个向量ALU指令。在Haswell架构上，其中5个uops（乘法和移位）只能在端口0上运行，因此无论如何展开这个算法，最多每5个周期就能实现一个向量（即每5/4个周期进行一次乘法）。
移位操作可以用
pshufb（端口5）来替代，以移动数据并在其中填充零。 （其他洗牌操作不支持用零替换而不是从输入中复制一个字节，并且输入中没有已知的零可以复制。）
paddq /
psubq可以在Haswell上的端口1/5或Skylake上的p015上运行。
Skylake在p01上运行pmuludq和立即计数向量移位，因此理论上可以每个max(5/2, 8/3, 11/4) = 11/4 = 2.75个周期处理一个向量。因此，它在总融合域uop吞吐量（包括2个向量加载和1个向量存储）上存在瓶颈。所以一点循环展开会有所帮助。可能由于不完美的调度而导致资源冲突，将其瓶颈限制在每个时钟周期略低于4个融合域uop。循环开销希望能在端口6上运行，该端口只能处理一些标量操作，包括add和compare-and-branch，将端口0/1/5留给向量ALU操作，因为它们接近饱和（8/3 = 2.666个时钟周期）。然而，加载/存储端口远未饱和。
所以，Skylake理论上可以在2.75个周期内处理一个向量（加上循环开销），或者在约0.7个周期内进行一次乘法，使用GCC自动向量化。相比之下，Haswell的最佳选择是理论上每1.2个周期进行一次标量操作，或者理论上每1.25个周期进行一次向量操作。然而，每1.2个周期的标量操作可能需要手动调优的汇编循环，因为编译器不知道如何使用单寄存器寻址模式进行存储，以及使用双寄存器寻址模式进行加载（dst + (src-dst) 并递增 dst）。
另外，如果你的数据不在L1缓存中，通过使用更少的指令完成任务可以让前端在执行单元之前提前进行加载，并在数据被需要之前开始加载。硬件预取不会跨越页面边界，因此对于大型数组，向量循环可能在实践中胜过标量循环，甚至对于较小的数组也可能如此。
AVX-512DQ引入了一个64bx64b->64b的向量乘法。
GCC可以使用它进行自动向量化，只需添加-mavx512dq选项。（实际上，使用-march=x86-64-v4，或者-march=native，或者-march=skylake-avx512等选项来启用其他功能并设置调优选项。）
.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（预计将成为Skylake多插槽Xeon（Purley）的一部分，大约在2017年）将通过更宽的向量提供比2倍更大的加速（如果这些指令每个时钟周期都能够进行流水线处理）。
更新：Skylake-AVX512（又称SKL-X或SKL-SP）以每1.5个周期运行VPMULLQ，适用于xmm、ymm或zmm向量。它是3个微操作，延迟为15个周期。（如果这不是AIDA结果中的测量故障，则zmm版本可能会额外增加1个周期的延迟。）
vpmullq比你用32位块构建的任何东西都要快，所以即使当前的CPU没有64位元素的向量乘法硬件，拥有这样的指令也是非常值得的。（可能它们使用FMA单元中的尾数乘法器。）
Zen 4将vpmullq作为单个微操作运行在两个端口中的任意一个，因此对于256位向量，每个时钟周期可以处理2个，对于512位向量，每个时钟周期可以处理1个。后来的英特尔CPU（如Alder Lake / Sapphire Rapids）仍然将其作为3个微操作运行在端口0/1上，因此每个时钟周期可以处理1.5个。https://uops.info/
__int64_t。只需#include <stdint.h>并使用int64_t，这样您的代码就具有可移植性了。幸运的是，64位乘法的低64位结果无论输入是有符号还是无符号都是相同的，因此该代码以任何方式都会给出相同的结果。 - Peter Cordes