你似乎在你的代码中假设long是64位，但也使用了__uint64_t。在32位中，x32 ABI，在Windows上，long是32位类型。你的标题提到了long long，但你的代码忽略了它。我想知道你的代码是否假设long是32位。
使用AVX 256加载，然后将指针别名化到__m256i上进行标量操作，这完全是在搬起石头砸自己的脚。gcc放弃了，给了你想要的糟糕代码：向量加载，然后是一堆extract和insert指令。你的写作方式意味着都vectors必须被解包，才能在标量中执行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）很好地自动向量化了一个简单的标量实现。当AVX 2不可用时，gcc仍然愚蠢地自动向量化了128 b个向量：在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 #,

型
如果你愿意的话，可以把它翻译回intrinsic，但是让编译器自动向量化会容易得多，我没有试图分析它，看看它是否是最优的。
如果你不经常使用-O3编译，你可以在循环之前使用#pragma omp simd（和-fopenmp）。
当然，与其使用标量结尾，不如对Gj_vec的最后32 B进行非对齐加载，并将其存储到Gi_vec的最后32 B中，这可能会与循环中的最后一次存储重叠。（如果数组小于32 B，则仍然需要标量回退。）
GCC 13和Clang 17 still both vectorize this way，使用另外两个vpmuludq而不是1x vpmulld来完成两个32位叉积，并稍微减少了 Shuffle 。但它们能够优化OP指针中的插入/提取-别名ptr_I[0] -= ptr_J[0] * q;，因此编译为与简单标量版本相同的asm循环。

AVX 2改进的vector intrinsic版本

从我对Z Boson的回答的评论中。基于Agner Fog的vector类库代码，在2023年沿着进行了各种改进，并修复了一个错误。（8年来，没有人注意到我将叉积添加到完整结果的底部，而不是添加到上半部分以匹配它们的位置值;感谢@Petr实际测试这个。）
Agner Fog的版本通过使用phadd + pshufd（我使用psllq/pand/paddd）保存了一条指令，但在shuffle端口上存在瓶颈。vpmulld将两个交叉产品打包在一起在Haswell和更高版本的Intel上是2个uop（与两个单独的vpmuludq指令相同，但输入的混洗较少），但在Zen 3和更高版本的AMD上，两个端口中的任何一个都只有1个uop。
由于其中一个操作数是常量，因此请确保将set1(q)作为b而不是a传递，以便可以提升“B_swap”shuffle。

// 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;  // same for signed or unsigned a,b since we aren't widening to 128
    // 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的吞吐量。用带有向量常量的vpshufb替换vpsllq（srli_epi64），或者用4替换vpslldq，并在添加后进行掩码，可以将端口0瓶颈降低到3个周期。）在Skylake和Rocket Lake等较新的CPU上，任何一种方式都可以使用，其中https://uica.uops.info/根据执行端口压力估计吞吐量瓶颈为2.67（只做乘法，没有循环或减法），有完美的调度。我用的是移位版本，所以它不需要额外的向量常数，并让shift/和并行运行，以获得更好的关键路径延迟，以防这是依赖链的一部分，不像在问题中，我们在加载和存储之间没有对每个元素做太多其他操作，因此无序的exec不必努力在迭代之间重叠。

我们可以通过执行sumcross = ((crossprod >> 32) + crossprod) << 32来完全避免向量常量。这仍然需要3条指令，但是AND变成了移位，所以它是另一个与乘法竞争的uop。（在Intel P核上）。它的关键路径延迟更差。但它避免了常量，使用更少的临时寄存器（出于同样的原因，它具有更差的ILP和更高的延迟），所以如果乘法作为一个更大的循环体的一部分，否则会耗尽寄存器，那么它可能是有用的。

AVX 1或SSE 4.1版本（每个向量两个元素）不会比64位标量imul r64, [mem]快。除非你已经有了向量数据，否则不要这样做，并希望得到一个向量。（这可能比提取到标量和回来更好。）或者在一个CPU上，有缓慢的imul r64, r64，也许像Silvermont-在这种情况下，_mm_add_epi32比_mm_add_epi64快。
add_epi64在任何AVX 2 CPU上都不会更慢或更大的代码大小，但是SSE版本在一些早期的CPU上。我尽可能使用add_epi32，因为它不会更慢，并且可能保存一些功率，或者在一些假设的未来CPU上可能具有更低的延迟。根据https://uops.info/，英特尔自Haswell和AMD自Zen 1运行vpaddd相同的vpaddq，包括桤木湖E核心。

GCC自动向量化代码的性能分析（非固有版本）

背景：参见Agner Fog's insn tables and microarch guide，以及x86标签wiki中的其他链接。
直到AVX 512（见下文），这可能只是略快于标量64位代码：imul r64, m64在Intel CPU上的吞吐量为每时钟一个（但在AMD Bulldozer-family上每4个时钟就有一个）. load/imul/sub with-memory-dest是英特尔CPU上的4个融合域uop（具有可以微融合的寻址模式，GCC未能使用该寻址模式）。流水线宽度是每时钟4个融合域uop，因此，即使是一个大展开也不能使其在每个时钟发出一次。有了足够的展开，我们将在加载/存储吞吐量上形成瓶颈。在Haswell上，每个时钟可以进行2次加载和一次存储，但是根据英特尔的手册，存储地址uop窃取加载端口将使吞吐量降低到大约81/96 = 84%。
所以也许Haswell最好的方法是加载并乘以标量，（2个uops），然后vmovq/pinsrq/vinserti128，这样你就可以用vpsubq做减法。这是8个uops来加载和乘以所有4个标量，7 shuffle uop，将数据放入__m256i（2（movq）+ 4（pinsrq是2 uop）+1 vinserti 128），还有3个uop用于执行vector load / vpsubq / vector store。（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的p015端口上运行。
Skylake在p01上运行pmuludq和立即计数向量移位，因此理论上它可以管理每个最大值一个向量的吞吐量（5/2，8/3，11/4）= 11/4 = 2.75个周期。因此，它会对总的融合域uop吞吐量造成瓶颈（包括2个向量加载和1个向量存储）。所以一点循环展开会有帮助。可能由于不完美的调度导致的资源冲突会使其瓶颈略低于4个融合-每个时钟的域uops。循环开销可以在端口6上运行，它只能处理一些标量操作，包括add和比较和分支，留下端口0/1/5用于向量ALU操作，因为它们接近饱和（8/3 = 2.666个时钟）。但是加载/存储端口远未接近饱和。
因此，Skylake理论上可以每2.75个周期管理一个向量（加上循环开销），或者使用GCC自动向量化每约0.7个循环执行一次乘法，而Haswell的最佳选择（理论上对于标量每1.2个周期一个，或者理论上对于向量每1.25个周期一个）。但是，标量每1.2个周期一个可能需要手动调优asm循环，因为编译器不知道如何使用一个寄存器的寻址模式来存储，以及两个寄存器的寻址模式来加载（dst + (src-dst)和增量dst）。

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

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

如果你添加-mavx512dq，GCC可以使用它自动矢量化。（在实践中，使用-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 #,

型
因此，AVX 512 DQ（预计将在2017年成为Skylake多插槽Xeon（Purley）的一部分）将给予比2倍大得多的加速（从更宽的向量），如果这些指令以每个时钟一个的速度流水线化。

更新：Skylake-AVX 512（也称为SKL-X或SKL-SP）每1.5个周期运行一次xmm，ymm或zmm向量的VPMULLQ。它是3个uop，15 c延迟。（如果不是测量故障in the AIDA results，zmm版本可能会有额外的1 c延迟。）
vpmullq比你用32位块构建的任何东西都要快得多，所以即使当前的CPU没有64位元素的向量乘法硬件，也非常值得为此编写一条指令。（假设它们在FMA单元中使用尾数乘数。）
Zen 4将vpmullq作为2个端口中的任一个的单个uop运行，因此256位向量的吞吐量为2/时钟，512位向量的吞吐量为1/时钟。后来的Intel CPU（如桤木Lake / Sapphire急流）仍然将其作为端口0/1的3个uop运行，因此吞吐量为1.5c。https://uops.info/

如果你对SIMD 64 bx 64 b到64 b（更低）操作感兴趣，这里有来自Agner Fog的Vector Class Library的AVX和AVX 2解决方案。我会用数组测试这些解决方案，看看它与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在循环中是常量，你不需要在每次迭代中重新计算一些东西，但你的编译器可能会发现这一点。