在 C 中使用 AVX 实现矩阵运算
Implementing matrix operation using AVX in C
我正在尝试使用 AVX 实现以下操作：
for (i=0; i<N; i++) {
for(j=0; j<N; j++) {
for (k=0; k<K; k++) {
d[i][j] += 2 * a[i][k] * ( b[k][j]- c[k]);
}
}
}
for (int i=0; i<N; i++){
f+= d[ind[i]][ind[i]]/2;
}
其中 d 是 NxN 矩阵，a 是 NxK，b 是 KxN，c 是长度为 K 的向量。它们都是双精度数。当然，所有数据都是对齐的，我正在使用 #pragma vector aligned 来帮助编译器 (gcc)。
我知道如何对一维数组使用 AVX 扩展，但对我来说用矩阵来做有点棘手。目前，我有以下内容，但没有得到正确的结果：
for (int i=0; i< floor (N/4); i++){
for (int j=0; j< floor (N/4); j++){
__m256d D, A, B, C;
D = _mm256_setzero_pd();
#pragma vector aligned
for (int k=0; k<K_MAX; k++){
A = _mm256_load_pd(a[i] + k*4);
B = _mm256_load_pd(b[k] + j*4);
C = _mm256_load_pd(c + 4*k);
B = _mm256_sub_pd(B, C);
A = _mm256_mul_pd(A, B);
D = _mm256_add_pd(_mm256_set1_pd(2.0), A);
_mm256_store_pd(d[i] + j*4, D);
}
}
}
for (int i=0; i<N; i++){
f+= d[ind[i]][ind[i]]/2;
}
希望有人能告诉我哪里错了
提前致谢。
NOTE: 我不太愿意介绍 OpenMP，只是使用 SIMD Intel 指令
假设 N 和 K 数字都相对较大（比硬件矢量大小 4 大很多），这是一种矢量化主循环的方法。未经测试。
主要思想是矢量化中间循环而不是内部循环。这样做有两个原因。
这样就避免了横向操作。当仅对内部循环进行向量化时，我们将不得不计算向量的水平和。
当加载 4 个连续的 k 值时，b[k][j] 加载具有不幸的 RAM 访问模式，需要 4 个单独的加载指令，或收集加载，这两种方法都相对较慢.为 4 个连续的 j 值加载元素是一个 full-vector 加载指令，非常有效，尤其是当您对齐输入时。
const int N_aligned = ( N / 4 ) * 4;
for( int i = 0; i < N; i++ )
{
int j = 0;
for( ; j < N_aligned; j += 4 )
{
// Load 4 scalars from d
__m256d dv = _mm256_loadu_pd( &d[ i ][ j ] );
// Run the inner loop which only loads from RAM but never stores any data
for( int k = 0; k < K; k++ )
{
__m256d av = _mm256_broadcast_sd( &a[ i ][ k ] );
__m256d bv = _mm256_loadu_pd( &b[ k ][ j ] );
__m256d cv = _mm256_broadcast_sd( &c[ k ] );
// dv += 2*av*( bv - cv )
__m256d t1 = _mm256_add_pd( av, av ); // 2*av
__m256d t2 = _mm256_sub_pd( bv, cv ); // bv - cv
dv = _mm256_fmadd_pd( t1, t2, dv );
}
// Store the updated 4 values
_mm256_storeu_pd( &d[ i ][ j ], dv );
}
// Handle remainder with scalar code
for( ; j < N; j++ )
{
double ds = d[ i ][ j ];
for( int k = 0; k < K; k++ )
ds += 2 * a[ i ][ k ] * ( b[ k ][ j ] - c[ k ] );
d[ i ][ j ] = ds;
}
}
如果你想进一步优化，尝试将内循环展开一个小因子，比如 2，使用 2 个用 _mm256_setzero_pd() 初始化的独立累加器，在循环之后添加它们。可能是在某些处理器上，此版本在 FMA 指令的延迟上停止，而不是使加载端口或 ALU 饱和。多个独立的累加器有时会有所帮助。
b[k][j] 是您的问题：元素 b[k + 0..3][j] 在内存中不连续。使用 SIMD（以合理/有用的方式）不是您可以放入经典的朴素 matmul 循环中的东西。请参阅 What Every Programmer Should Know About Memory? - 有一个附录，其中包含一个 SSE2 matmul（带有 cache-blocking）的示例，它显示了如何以不同的顺序执行操作，即 SIMD-friendly。
展示了如何通过对中间循环 j 进行矢量化来进行矢量化。但这留下了一个相对较差的内存访问模式，并且循环内有 3 个加载 + 3 个 ALU 操作。 （这个答案开始是对它的评论，请参阅它以获取我正在谈论的代码并提出更改建议。）
循环反转应该可以像 inner-most 循环 那样执行 j。这意味着在 inner-most 循环中为 d[i][j] += ... 做存储，但 OTOH 它在 2 * a[i][k] * ( b[k][j]- c[k] ) 中产生更多循环不变量，因此您可以有效地转换为 d[i][j] += (2*a_ik) * b[k][j] - (2*a_ik*c_k)，即一个 VFMSUBPD 和一个每个加载和存储的 VADDPD。 （bv 加载折叠到 FMSUB 作为内存源操作数，dv 加载折叠到 VADDPD，所以希望 front-end 只有 3 微指令，包括一个单独的存储，而不是包括循环开销。）
编译器必须在 Intel CPU 的端口 7 上展开 avoid an indexed addressing mode so the store-address uop can stay micro-fused 和 运行（Haswell 到 Skylake-family），而不是与两个负载竞争。 Ice Lake 没有这个问题，有两个完全独立的 store-AGUs 与两个负载 AGU 分开。但可能仍需要展开一些循环以避免 front-end 瓶颈。
这是一个例子，未经测试（Soonts 提供的原始版本，谢谢）。它以不同的方式优化了循环中的 2 FP 数学运算：简单地将 2*a 提升到循环之外，对 dv += (2av)*(sub_result) 执行 SUB 然后 FMA。但是 bv 不能作为 vsubpd 的源操作数，因为我们需要 bv - cv。但是我们可以通过否定 cv 以允许在内循环中使用 (-cv) + bv 来解决这个问题，并使用 bv 作为内存源操作数。有时编译器会为你做这样的事情，但在这里他们似乎没有，所以我手动做了。否则我们会得到一个单独的 vmovupd 负载通过 front-end.
#include <stdint.h>
#include <stdlib.h>
#include <immintrin.h>
// This double [N][N] C99 VLA syntax isn't portable to C++ even with GNU extensions
// restrict tells the compiler the output doesn't overlap with any of the inputs
void matop(size_t N, size_t K, double d[restrict N][N], const double a[restrict N][K], const double b[restrict K][N], const double c[restrict K])
{
for( size_t i = 0; i < N; i++ ) {
// loop-invariant pointers for this outer iteration
//double* restrict rowDi = &d[ i ][ 0 ];
const double* restrict rowAi = &a[ i ][ 0 ];
for( size_t k = 0; k < K; k++ ) {
const double* restrict rowBk = &b[ k ][ 0 ];
double* restrict rowDi = &d[ i ][ 0 ];
#if 0 // pure scalar
// auto-vectorizes ok; still a lot of extra checking outside outermost loop even with restrict
for (size_t j=0 ; j<N ; j++){
rowDi[j] += 2*rowAi[k] * (rowBk[j] - c[k]);
}
#else // SIMD inner loop with cleanup
// *** TODO: unroll over 2 or 3 i values
// and maybe also 2 or 3 k values, to reuse each bv a few times while it's loaded.
__m256d av = _mm256_broadcast_sd( rowAi + k );
av = _mm256_add_pd( av, av ); // 2*a[ i ][ k ] broadcasted
const __m256d cv = _mm256_broadcast_sd( &c[ k ] );
const __m256d minus_ck = _mm256_xor_pd(cv, _mm256_set1_pd(-0.0)); // broadcasted -c[k]
//const size_t N_aligned = ( (size_t)N / 4 ) * 4;
size_t N_aligned = N & -4; // round down to a multiple of 4 j iterations
const double* endBk = rowBk + N_aligned;
//for( ; j < N_aligned; j += 4 )
for ( ; rowBk != endBk ; rowBk += 4, rowDi += 4) { // coax GCC into using pointer-increments in the asm, instead of j+=4
// Load the output vector to update
__m256d dv = _mm256_loadu_pd( rowDi );
// Update with FMA
__m256d bv = _mm256_loadu_pd( rowBk );
__m256d t2 = _mm256_add_pd( minus_ck, bv ); // bv - cv
dv = _mm256_fmadd_pd( av, t2, dv );
// Store back to the same address
_mm256_storeu_pd( rowDi, dv );
}
// rowDi and rowBk point to the double after the last full vector
// The remainder, if you can't pad your rows to a multiple of 4 and step on that padding
for(int j=0 ; j < (N&3); j++ )
rowDi[ j ] += _mm256_cvtsd_f64( av ) * ( rowBk[ j ] + _mm256_cvtsd_f64( minus_ck ) );
#endif
}
}
}
没有展开（https://godbolt.org/z/6WeYKbnYY），GCC11 的内部循环 asm 看起来像这样，所有 single-uop 指令可以保留 micro-fused 即使在 back-end on Haswell 和更高版本中.
.L7: # do{
vaddpd ymm0, ymm2, YMMWORD PTR [rax] # -c[k] + rowBk[0..3]
add rax, 32 # rowBk += 4
add rdx, 32 # rowDi += 4
vfmadd213pd ymm0, ymm1, YMMWORD PTR [rdx-32] # fma(2aik, Bkj-ck, Dij)
vmovupd YMMWORD PTR [rdx-32], ymm0 # store FMA result
cmp rcx, rax
jne .L7 # }while(p != endp)
但它总共有 6 个微指令，其中 3 个循环开销（指针递增和融合的 cmp+jne），因此通过 Skylake 的 Haswell 只能 运行 它每 1.5 个时钟迭代 1 次，瓶颈在 4- front-end 的广泛发行阶段。 （这不会让 OoO exec 在执行指针增量和循环分支方面取得领先，到 而 back-end 仍在咀嚼较旧的负载和 FP 数学。）
所以循环展开应该会有所帮助，因为我们设法诱使 GCC 使用索引寻址模式。如果没有它，它在 Intel Haswell/Skylake CPU 上使用 AVX 代码相对没用，每个 vaddpd ymm5, ymm4, [rax + r14] 解码为 1 micro-fused uop，但在 back-end 中分解成 2 个，无济于事我们通过 front-end 的最窄部分获得更多工作。 （很像如果我们使用单独的 vmovupd 加载，就像我们使用 _mm256_sub_pd(bv, cv) 而不是 add(bv, -cv) 一样。）
vmovupd ymmword ptr [rbp + r14], ymm5 存储保持 micro-fused 但不能 运行 在端口 7 上，限制我们每个时钟总共 2 次内存操作（最多 1 次可以一个商店。）所以最好的情况是每个向量 1.5 个周期。
使用 GCC 和 clang -O3 -march=skylake -funroll-loops 在 https://godbolt.org/z/rd3rn9zor 上编译。 GCC 实际上使用指针增量，负载折叠成 8x vaddpd 和 8x vfmadd213pd。但是 clang 使用索引寻址模式并且不会展开。 （您可能不希望整个程序使用 -funroll-loops，因此要么单独编译它，要么手动展开。GCC 的展开完全剥离了一个在进入实际 SIMD 循环之前进行 0..7 次向量迭代的序言，所以它相当咄咄逼人。）
GCC 的 loop-unrolling 在这里对于大 N 看起来很有用，分摊指针增量和多个向量上的循环开销。 （例如，GCC 不知道如何在点积中为 FP dep 链发明多个累加器，这使得它的展开在那种情况下毫无用处，这与 clang 不同。）
不幸的是，clang 没有为我们展开内部循环，但它确实以一种有趣的方式使用 vmaskmovpd 进行清理。
我们 使用单独的循环计数器进行清理可能会更好，这样编译器可以轻松地证明 trip-count 用于清理的是 0..3, 所以它不会尝试 auto-vectorize 与 YMM.
另一种方法，使用实际的 j 变量进行内部循环及其清理，更像是 Soonts 的编辑。 IIRC，编译器 did 尝试 auto-vectorize 清理，浪费代码大小和一些 always-false 分支。
size_t j = 0; // used for cleanup loop after
for( ; j < N_aligned; j += 4 )
{
// Load the output vector to update
__m256d dv = _mm256_loadu_pd( rowDi + j );
// Update with FMA
__m256d bv = _mm256_loadu_pd( rowBk + j );
__m256d t2 = _mm256_sub_pd( bv, cv ); // bv - cv
dv = _mm256_fmadd_pd( av, t2, dv );
// Store back to the same address
_mm256_storeu_pd( rowDi + j, dv );
}
// The remainder, if you can't pad your rows to a multiple of 4
for( ; j < N; j++ )
rowDi[ j ] += _mm256_cvtsd_f64( av ) * ( rowBk[ j ] - _mm256_cvtsd_f64( cv ) );
现代 CPU (https://agner.org/optimize/ and https://uops.info/) 的加载和存储与 FP 数学的混合相当好，尤其是英特尔，我们可以进行 2 次加载 和 1 次存储。我认为 Zen 2 或 3 也可以做 2 loads + 1 store。它需要命中 L1d 缓存以维持这种吞吐量， 尽管。 （即便如此，英特尔的优化手册说 Skylake 上的最大持续 L1d 带宽小于所需的全部 96 bytes/cycle。更像是 80 年代中期的 IIRC，所以我们不能完全期望每个周期一个结果向量，即使有足够的展开来避免 front-end 瓶颈。）
没有延迟瓶颈，因为我们在每次迭代中都转向一个新的 dv，而不是在循环迭代中累积任何东西。
这样做的另一个好处是 对 d[i][j] 和 b[k][j] 的内存访问是顺序的，[=192] 中没有其他内存访问=] 循环。 （中间循环将执行 a[i][k] 和 c[k] 的 broadcast-loads。如果内部循环驱逐太多，这些似乎可能 cache-miss；随着外部循环的展开，一个SIMD 加载和一些改组可能会有所帮助，但可能 cache-blocking 会避免这种需要。）
针对不同的 b[k] 行重复遍历相同的 d[i] 行为我们提供了我们正在修改的部分的局部性（即使用 k 作为中间循环，保持 i 作为 outer-most.) 使用 k 作为外层循环，我们将在整个 d[0..N-1][0..N-1] 中循环 K 次，可能需要写 + 读每一个都传递到任何级别的缓存或内存可以容纳它。
但如果每一行真的很长，你仍然想要 cache-block，这样你就可以避免缓存未命中带来所有 b[][]从 DRAM 输入 N 次。并避免驱逐你接下来要 broadcast-load 的东西。
更智能的展开：迈向 cache-blocking
的第一步
如果我们对每个数据向量做更多的处理，那么上述一些与最大化 load/store 执行单元吞吐量和要求编译器使用 non-indexed 寻址模式有关的问题可能会消失已加载。
例如，我们可以处理 2、3 或 4，而不是只处理 d[][] 的一行。然后每个 (rowBk[j] - c[k]) 结果都可以使用那么多次（用d[i+unroll][j + 0..vec] 向量的不同 2aik)。
我们还可以加载几个不同的 (rowBk+K*0..unroll)[j+0..3]，每个都有相应的 minus_ck0、minus_ck1 等（或者保留一个向量数组；只要它很小并且编译器有足够的寄存器，元素将不存在于内存中。）
同时在寄存器中使用多个 bv-cv 和 dv 向量，我们可以在不增加 FP 工作总量的情况下为每次加载执行更多的 FMA。但是，它需要更多的常量寄存器，否则我们可能会通过强制重新加载来破坏目的。
d[i][j] += (2*a_ik) * b[k][j] - (2*a_ik*c_k) 转换在这里没有用；我们希望将 bv-cv 与 i 分开，以便我们可以将该结果重新用作不同 FMA 的输入。
b[k][j]+(-c[k]) 仍然可以从带有 vaddpd 的负载的 micro-fusion 中获益，因此理想情况下它仍会使用指针增量，但 front-end 可能不会瓶颈了。
这个不要做得太过分了；太多的内存输入流可能是缓存冲突未命中的问题，尤其是对于可能产生别名的某些 N 值，以及硬件预取跟踪它们。 （尽管据说英特尔的 L2 流媒体每 4k 页跟踪 1 个前向流和 1 个后向流，IIRC。）可能大约 4 到 8 个 ish 流是可以的。但是如果 d[][] 在 L1d 中没有丢失，那么它就不是真正的内存输入流。但是，您不希望 b[][] 输入行逐出 d 数据，因为您将重复循环 2 到 4 行 d 数据。
相比之下：Soonts 循环 - 清理频率较低，但内存访问模式较差。
Soonts 的当前循环有 3 个负载和 3 个 ALU 操作并不理想，尽管如果它们命中缓存，每个 FMA 操作 1 个负载已经可以了（大多数现代 CPU 每个时钟可以做 2 个，尽管 AMD Zen 也可以做 2 FP 与 mul/fma 并行添加）。如果额外的 ALU 操作是一个瓶颈，我们可以 pre-multiply a[][] 乘以 2 一次，只需要 O(N*K) 工作与 O(N^2*K) 即时完成。但这可能不是瓶颈，因此不值得。
更重要的是，Soonts 当前答案中的内存访问模式是一次向前循环 1 次，用于 c[k] 和 a[i][k] 的广播负载，这很好，但是 bv = _mm256_loadu_pd of b[k][j + 0..3] 不幸跨栏。
如果你要按照 Soonts 的建议展开，不要只为一个 dv 做两个 dep 链，至少做两个向量，d[i][j + 0..3] 和 4..7 所以您使用每个 b[k][j] 触摸的整个 64 字节（完整缓存行）。或一对 cache-lines 的四个向量。 （Intel CPU 至少使用 adjacent-line 预取器，它喜欢完成 128 字节的 aligned 对缓存行，所以你会受益于 b[][] 128。或至少 64，并从 adjacent-line 预取中获得一些好处。
如果 b[][] 的垂直切片适合某个级别的缓存（连同您当前正在累积的 d[i][] 的行），则下一个 gr 的下一步最多的列可以从该预取和位置中受益。如果没有，充分利用你接触的线更重要，这样他们以后就不必再被拉进来了。
因此，使用 Soonts 的矢量化策略，对于 L1d 缓存不适合的大问题，确保 b 的行按 64 位对齐可能很好，即使这意味着在末尾进行填充每一行。 （存储几何不必匹配实际的矩阵维度；您分别传递 N 和 row_stride。您将一个用于索引计算，另一个用于循环边界。）
我正在尝试使用 AVX 实现以下操作：
for (i=0; i<N; i++) {
for(j=0; j<N; j++) {
for (k=0; k<K; k++) {
d[i][j] += 2 * a[i][k] * ( b[k][j]- c[k]);
}
}
}
for (int i=0; i<N; i++){
f+= d[ind[i]][ind[i]]/2;
}
其中 d 是 NxN 矩阵，a 是 NxK，b 是 KxN，c 是长度为 K 的向量。它们都是双精度数。当然，所有数据都是对齐的，我正在使用 #pragma vector aligned 来帮助编译器 (gcc)。
我知道如何对一维数组使用 AVX 扩展，但对我来说用矩阵来做有点棘手。目前，我有以下内容，但没有得到正确的结果：
for (int i=0; i< floor (N/4); i++){
for (int j=0; j< floor (N/4); j++){
__m256d D, A, B, C;
D = _mm256_setzero_pd();
#pragma vector aligned
for (int k=0; k<K_MAX; k++){
A = _mm256_load_pd(a[i] + k*4);
B = _mm256_load_pd(b[k] + j*4);
C = _mm256_load_pd(c + 4*k);
B = _mm256_sub_pd(B, C);
A = _mm256_mul_pd(A, B);
D = _mm256_add_pd(_mm256_set1_pd(2.0), A);
_mm256_store_pd(d[i] + j*4, D);
}
}
}
for (int i=0; i<N; i++){
f+= d[ind[i]][ind[i]]/2;
}
希望有人能告诉我哪里错了
提前致谢。
NOTE: 我不太愿意介绍 OpenMP，只是使用 SIMD Intel 指令
假设 N 和 K 数字都相对较大（比硬件矢量大小 4 大很多），这是一种矢量化主循环的方法。未经测试。
主要思想是矢量化中间循环而不是内部循环。这样做有两个原因。
这样就避免了横向操作。当仅对内部循环进行向量化时，我们将不得不计算向量的水平和。
当加载 4 个连续的
k值时，b[k][j]加载具有不幸的 RAM 访问模式，需要 4 个单独的加载指令，或收集加载，这两种方法都相对较慢.为 4 个连续的j值加载元素是一个 full-vector 加载指令，非常有效，尤其是当您对齐输入时。
const int N_aligned = ( N / 4 ) * 4;
for( int i = 0; i < N; i++ )
{
int j = 0;
for( ; j < N_aligned; j += 4 )
{
// Load 4 scalars from d
__m256d dv = _mm256_loadu_pd( &d[ i ][ j ] );
// Run the inner loop which only loads from RAM but never stores any data
for( int k = 0; k < K; k++ )
{
__m256d av = _mm256_broadcast_sd( &a[ i ][ k ] );
__m256d bv = _mm256_loadu_pd( &b[ k ][ j ] );
__m256d cv = _mm256_broadcast_sd( &c[ k ] );
// dv += 2*av*( bv - cv )
__m256d t1 = _mm256_add_pd( av, av ); // 2*av
__m256d t2 = _mm256_sub_pd( bv, cv ); // bv - cv
dv = _mm256_fmadd_pd( t1, t2, dv );
}
// Store the updated 4 values
_mm256_storeu_pd( &d[ i ][ j ], dv );
}
// Handle remainder with scalar code
for( ; j < N; j++ )
{
double ds = d[ i ][ j ];
for( int k = 0; k < K; k++ )
ds += 2 * a[ i ][ k ] * ( b[ k ][ j ] - c[ k ] );
d[ i ][ j ] = ds;
}
}
如果你想进一步优化，尝试将内循环展开一个小因子，比如 2，使用 2 个用 _mm256_setzero_pd() 初始化的独立累加器，在循环之后添加它们。可能是在某些处理器上，此版本在 FMA 指令的延迟上停止，而不是使加载端口或 ALU 饱和。多个独立的累加器有时会有所帮助。
b[k][j] 是您的问题：元素 b[k + 0..3][j] 在内存中不连续。使用 SIMD（以合理/有用的方式）不是您可以放入经典的朴素 matmul 循环中的东西。请参阅 What Every Programmer Should Know About Memory? - 有一个附录，其中包含一个 SSE2 matmul（带有 cache-blocking）的示例，它显示了如何以不同的顺序执行操作，即 SIMD-friendly。
j 进行矢量化来进行矢量化。但这留下了一个相对较差的内存访问模式，并且循环内有 3 个加载 + 3 个 ALU 操作。 （这个答案开始是对它的评论，请参阅它以获取我正在谈论的代码并提出更改建议。）
循环反转应该可以像 inner-most 循环 那样执行 j。这意味着在 inner-most 循环中为 d[i][j] += ... 做存储，但 OTOH 它在 2 * a[i][k] * ( b[k][j]- c[k] ) 中产生更多循环不变量，因此您可以有效地转换为 d[i][j] += (2*a_ik) * b[k][j] - (2*a_ik*c_k)，即一个 VFMSUBPD 和一个每个加载和存储的 VADDPD。 （bv 加载折叠到 FMSUB 作为内存源操作数，dv 加载折叠到 VADDPD，所以希望 front-end 只有 3 微指令，包括一个单独的存储，而不是包括循环开销。）
编译器必须在 Intel CPU 的端口 7 上展开 avoid an indexed addressing mode so the store-address uop can stay micro-fused 和 运行（Haswell 到 Skylake-family），而不是与两个负载竞争。 Ice Lake 没有这个问题，有两个完全独立的 store-AGUs 与两个负载 AGU 分开。但可能仍需要展开一些循环以避免 front-end 瓶颈。
这是一个例子，未经测试（Soonts 提供的原始版本，谢谢）。它以不同的方式优化了循环中的 2 FP 数学运算：简单地将 2*a 提升到循环之外，对 dv += (2av)*(sub_result) 执行 SUB 然后 FMA。但是 bv 不能作为 vsubpd 的源操作数，因为我们需要 bv - cv。但是我们可以通过否定 cv 以允许在内循环中使用 (-cv) + bv 来解决这个问题，并使用 bv 作为内存源操作数。有时编译器会为你做这样的事情，但在这里他们似乎没有，所以我手动做了。否则我们会得到一个单独的 vmovupd 负载通过 front-end.
#include <stdint.h>
#include <stdlib.h>
#include <immintrin.h>
// This double [N][N] C99 VLA syntax isn't portable to C++ even with GNU extensions
// restrict tells the compiler the output doesn't overlap with any of the inputs
void matop(size_t N, size_t K, double d[restrict N][N], const double a[restrict N][K], const double b[restrict K][N], const double c[restrict K])
{
for( size_t i = 0; i < N; i++ ) {
// loop-invariant pointers for this outer iteration
//double* restrict rowDi = &d[ i ][ 0 ];
const double* restrict rowAi = &a[ i ][ 0 ];
for( size_t k = 0; k < K; k++ ) {
const double* restrict rowBk = &b[ k ][ 0 ];
double* restrict rowDi = &d[ i ][ 0 ];
#if 0 // pure scalar
// auto-vectorizes ok; still a lot of extra checking outside outermost loop even with restrict
for (size_t j=0 ; j<N ; j++){
rowDi[j] += 2*rowAi[k] * (rowBk[j] - c[k]);
}
#else // SIMD inner loop with cleanup
// *** TODO: unroll over 2 or 3 i values
// and maybe also 2 or 3 k values, to reuse each bv a few times while it's loaded.
__m256d av = _mm256_broadcast_sd( rowAi + k );
av = _mm256_add_pd( av, av ); // 2*a[ i ][ k ] broadcasted
const __m256d cv = _mm256_broadcast_sd( &c[ k ] );
const __m256d minus_ck = _mm256_xor_pd(cv, _mm256_set1_pd(-0.0)); // broadcasted -c[k]
//const size_t N_aligned = ( (size_t)N / 4 ) * 4;
size_t N_aligned = N & -4; // round down to a multiple of 4 j iterations
const double* endBk = rowBk + N_aligned;
//for( ; j < N_aligned; j += 4 )
for ( ; rowBk != endBk ; rowBk += 4, rowDi += 4) { // coax GCC into using pointer-increments in the asm, instead of j+=4
// Load the output vector to update
__m256d dv = _mm256_loadu_pd( rowDi );
// Update with FMA
__m256d bv = _mm256_loadu_pd( rowBk );
__m256d t2 = _mm256_add_pd( minus_ck, bv ); // bv - cv
dv = _mm256_fmadd_pd( av, t2, dv );
// Store back to the same address
_mm256_storeu_pd( rowDi, dv );
}
// rowDi and rowBk point to the double after the last full vector
// The remainder, if you can't pad your rows to a multiple of 4 and step on that padding
for(int j=0 ; j < (N&3); j++ )
rowDi[ j ] += _mm256_cvtsd_f64( av ) * ( rowBk[ j ] + _mm256_cvtsd_f64( minus_ck ) );
#endif
}
}
}
没有展开（https://godbolt.org/z/6WeYKbnYY），GCC11 的内部循环 asm 看起来像这样，所有 single-uop 指令可以保留 micro-fused 即使在 back-end on Haswell 和更高版本中.
.L7: # do{
vaddpd ymm0, ymm2, YMMWORD PTR [rax] # -c[k] + rowBk[0..3]
add rax, 32 # rowBk += 4
add rdx, 32 # rowDi += 4
vfmadd213pd ymm0, ymm1, YMMWORD PTR [rdx-32] # fma(2aik, Bkj-ck, Dij)
vmovupd YMMWORD PTR [rdx-32], ymm0 # store FMA result
cmp rcx, rax
jne .L7 # }while(p != endp)
但它总共有 6 个微指令，其中 3 个循环开销（指针递增和融合的 cmp+jne），因此通过 Skylake 的 Haswell 只能 运行 它每 1.5 个时钟迭代 1 次，瓶颈在 4- front-end 的广泛发行阶段。 （这不会让 OoO exec 在执行指针增量和循环分支方面取得领先，到
所以循环展开应该会有所帮助，因为我们设法诱使 GCC 使用索引寻址模式。如果没有它，它在 Intel Haswell/Skylake CPU 上使用 AVX 代码相对没用，每个 vaddpd ymm5, ymm4, [rax + r14] 解码为 1 micro-fused uop，但在 back-end 中分解成 2 个，无济于事我们通过 front-end 的最窄部分获得更多工作。 （很像如果我们使用单独的 vmovupd 加载，就像我们使用 _mm256_sub_pd(bv, cv) 而不是 add(bv, -cv) 一样。）
vmovupd ymmword ptr [rbp + r14], ymm5 存储保持 micro-fused 但不能 运行 在端口 7 上，限制我们每个时钟总共 2 次内存操作（最多 1 次可以一个商店。）所以最好的情况是每个向量 1.5 个周期。
使用 GCC 和 clang -O3 -march=skylake -funroll-loops 在 https://godbolt.org/z/rd3rn9zor 上编译。 GCC 实际上使用指针增量，负载折叠成 8x vaddpd 和 8x vfmadd213pd。但是 clang 使用索引寻址模式并且不会展开。 （您可能不希望整个程序使用 -funroll-loops，因此要么单独编译它，要么手动展开。GCC 的展开完全剥离了一个在进入实际 SIMD 循环之前进行 0..7 次向量迭代的序言，所以它相当咄咄逼人。）
GCC 的 loop-unrolling 在这里对于大 N 看起来很有用，分摊指针增量和多个向量上的循环开销。 （例如，GCC 不知道如何在点积中为 FP dep 链发明多个累加器，这使得它的展开在那种情况下毫无用处，这与 clang 不同。）
不幸的是，clang 没有为我们展开内部循环，但它确实以一种有趣的方式使用 vmaskmovpd 进行清理。
我们 使用单独的循环计数器进行清理可能会更好，这样编译器可以轻松地证明 trip-count 用于清理的是 0..3, 所以它不会尝试 auto-vectorize 与 YMM.
另一种方法，使用实际的 j 变量进行内部循环及其清理，更像是 Soonts 的编辑。 IIRC，编译器 did 尝试 auto-vectorize 清理，浪费代码大小和一些 always-false 分支。
size_t j = 0; // used for cleanup loop after
for( ; j < N_aligned; j += 4 )
{
// Load the output vector to update
__m256d dv = _mm256_loadu_pd( rowDi + j );
// Update with FMA
__m256d bv = _mm256_loadu_pd( rowBk + j );
__m256d t2 = _mm256_sub_pd( bv, cv ); // bv - cv
dv = _mm256_fmadd_pd( av, t2, dv );
// Store back to the same address
_mm256_storeu_pd( rowDi + j, dv );
}
// The remainder, if you can't pad your rows to a multiple of 4
for( ; j < N; j++ )
rowDi[ j ] += _mm256_cvtsd_f64( av ) * ( rowBk[ j ] - _mm256_cvtsd_f64( cv ) );
现代 CPU (https://agner.org/optimize/ and https://uops.info/) 的加载和存储与 FP 数学的混合相当好，尤其是英特尔，我们可以进行 2 次加载 和 1 次存储。我认为 Zen 2 或 3 也可以做 2 loads + 1 store。它需要命中 L1d 缓存以维持这种吞吐量， 尽管。 （即便如此，英特尔的优化手册说 Skylake 上的最大持续 L1d 带宽小于所需的全部 96 bytes/cycle。更像是 80 年代中期的 IIRC，所以我们不能完全期望每个周期一个结果向量，即使有足够的展开来避免 front-end 瓶颈。）
没有延迟瓶颈，因为我们在每次迭代中都转向一个新的 dv，而不是在循环迭代中累积任何东西。
这样做的另一个好处是 对 d[i][j] 和 b[k][j] 的内存访问是顺序的，[=192] 中没有其他内存访问=] 循环。 （中间循环将执行 a[i][k] 和 c[k] 的 broadcast-loads。如果内部循环驱逐太多，这些似乎可能 cache-miss；随着外部循环的展开，一个SIMD 加载和一些改组可能会有所帮助，但可能 cache-blocking 会避免这种需要。）
针对不同的 b[k] 行重复遍历相同的 d[i] 行为我们提供了我们正在修改的部分的局部性（即使用 k 作为中间循环，保持 i 作为 outer-most.) 使用 k 作为外层循环，我们将在整个 d[0..N-1][0..N-1] 中循环 K 次，可能需要写 + 读每一个都传递到任何级别的缓存或内存可以容纳它。
但如果每一行真的很长，你仍然想要 cache-block，这样你就可以避免缓存未命中带来所有 b[][]从 DRAM 输入 N 次。并避免驱逐你接下来要 broadcast-load 的东西。
更智能的展开：迈向 cache-blocking
的第一步如果我们对每个数据向量做更多的处理，那么上述一些与最大化 load/store 执行单元吞吐量和要求编译器使用 non-indexed 寻址模式有关的问题可能会消失已加载。
例如，我们可以处理 2、3 或 4，而不是只处理 d[][] 的一行。然后每个 (rowBk[j] - c[k]) 结果都可以使用那么多次（用d[i+unroll][j + 0..vec] 向量的不同 2aik)。
我们还可以加载几个不同的 (rowBk+K*0..unroll)[j+0..3]，每个都有相应的 minus_ck0、minus_ck1 等（或者保留一个向量数组；只要它很小并且编译器有足够的寄存器，元素将不存在于内存中。）
同时在寄存器中使用多个 bv-cv 和 dv 向量，我们可以在不增加 FP 工作总量的情况下为每次加载执行更多的 FMA。但是，它需要更多的常量寄存器，否则我们可能会通过强制重新加载来破坏目的。
d[i][j] += (2*a_ik) * b[k][j] - (2*a_ik*c_k) 转换在这里没有用；我们希望将 bv-cv 与 i 分开，以便我们可以将该结果重新用作不同 FMA 的输入。
b[k][j]+(-c[k]) 仍然可以从带有 vaddpd 的负载的 micro-fusion 中获益，因此理想情况下它仍会使用指针增量，但 front-end 可能不会瓶颈了。
这个不要做得太过分了；太多的内存输入流可能是缓存冲突未命中的问题，尤其是对于可能产生别名的某些 N 值，以及硬件预取跟踪它们。 （尽管据说英特尔的 L2 流媒体每 4k 页跟踪 1 个前向流和 1 个后向流，IIRC。）可能大约 4 到 8 个 ish 流是可以的。但是如果 d[][] 在 L1d 中没有丢失，那么它就不是真正的内存输入流。但是，您不希望 b[][] 输入行逐出 d 数据，因为您将重复循环 2 到 4 行 d 数据。
相比之下：Soonts 循环 - 清理频率较低，但内存访问模式较差。
Soonts 的当前循环有 3 个负载和 3 个 ALU 操作并不理想，尽管如果它们命中缓存，每个 FMA 操作 1 个负载已经可以了（大多数现代 CPU 每个时钟可以做 2 个，尽管 AMD Zen 也可以做 2 FP 与 mul/fma 并行添加）。如果额外的 ALU 操作是一个瓶颈，我们可以 pre-multiply a[][] 乘以 2 一次，只需要 O(N*K) 工作与 O(N^2*K) 即时完成。但这可能不是瓶颈，因此不值得。
更重要的是，Soonts 当前答案中的内存访问模式是一次向前循环 1 次，用于 c[k] 和 a[i][k] 的广播负载，这很好，但是 bv = _mm256_loadu_pd of b[k][j + 0..3] 不幸跨栏。
如果你要按照 Soonts 的建议展开，不要只为一个 dv 做两个 dep 链，至少做两个向量，d[i][j + 0..3] 和 4..7 所以您使用每个 b[k][j] 触摸的整个 64 字节（完整缓存行）。或一对 cache-lines 的四个向量。 （Intel CPU 至少使用 adjacent-line 预取器，它喜欢完成 128 字节的 aligned 对缓存行，所以你会受益于 b[][] 128。或至少 64，并从 adjacent-line 预取中获得一些好处。
如果 b[][] 的垂直切片适合某个级别的缓存（连同您当前正在累积的 d[i][] 的行），则下一个 gr 的下一步最多的列可以从该预取和位置中受益。如果没有，充分利用你接触的线更重要，这样他们以后就不必再被拉进来了。
因此，使用 Soonts 的矢量化策略，对于 L1d 缓存不适合的大问题，确保 b 的行按 64 位对齐可能很好，即使这意味着在末尾进行填充每一行。 （存储几何不必匹配实际的矩阵维度；您分别传递 N 和 row_stride。您将一个用于索引计算，另一个用于循环边界。）