相关论文是here。我试图重现Kazushige Goto的开创性论文，通过将其衰减到gepp（通用面板 - 面板）和gebp（通用块面板）乘法的子程序来快速矩阵乘法，这显然是gemm的最快构建块。我编写了下面的代码来测试它并使用-O3标志，我看到我的代码的性能实际上比天真矩阵乘法更差：
(~0.5x increase)
Time elapsed : 3.82941 <-- naive
Time elapsed : 6.21072 <-- more complex gebp subroutine
然而，没有-O3标志，我们看到速度确实比天真版本更快：
(~4x increase)
Time elapsed : 53.4537 <-- naive
Time elapsed : 15.603 <-- more complex gebp subroutine
根据@ ztik的建议，我在没有-mavx2 -O3标志的情况下尝试了它，并添加了-O2，它显示了类似的结果，没有任何优化标志：
(~4x increase)
Time elapsed : 26.4217 <-- naive
Time elapsed : 6.42583 <-- more complex gebp subroutine
我很清楚-O3在gcc中启用了数百个优化标志，每个标志都会比复杂的Goto纸张变体更多地单独提高朴素方法的性能。然而，MKL，ATLAS等是使用Goto方法的一些着名的BLAS变体（并且比天真更快的数量级），尽管使用汇编内核而不是C ++代码。
这是预期的，如果是这样，我怎样才能获得性能提升（除了必须编写自展开的程序集）？我是否只是编写了一些可怕的代码，其中包含某些不会滥用缓存的错误？
我正在运行带有AVX2指令的macbook pro，带有intel生成：
Intel(R) Core(TM) i7-6660U CPU @ 2.40GHz
我使用g++-8编译下面的代码，它目前是v8.2.0。
可重现的代码如下：
Makefile：
CC=g++-8
FLAGS=-std=c++17 -ffast-math -mavx2 -O3
run: main
./main
main: main.cpp
$(CC) -o main main.cpp $(FLAGS)
main.cpp，包含一些用于基准测试的实用程序：
#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>
#define TIME_IT(EXPR) \
{ \
std::clock_t __start; \
double __duration; \
__start = std::clock(); \
EXPR; \
__duration = (std::clock() - __start) / (double) CLOCKS_PER_SEC; \
std::cout << "Time elapsed : " << __duration << std::endl; \
}
static constexpr auto X_SIZE = 2048;
static constexpr auto Y_SIZE = 1024;
static constexpr auto Z_SIZE = 2048;
// = 32 floats
static constexpr auto BLK_BYTES = 128;
static constexpr auto BLK_SIZE = BLK_BYTES / 4;
template <size_t row, size_t mid, size_t col>
void initialize_matrices(float (&a)[row][mid], float (&b)[mid][col], float (&c)[row][col]){
// Initialize matrices
for(auto i = 0; i < row; i++){
for(auto j = 0; j < mid; j++){
a[i][j] = ((float) (i*Y_SIZE + j)) / (row * mid);
}
}
for(auto i = 0; i < mid; i++){
for(auto j = 0; j < col; j++){
b[i][j] = ((float) (i*Z_SIZE + j)) / (mid * col);
}
}
for(auto i = 0; i < row; i++){
for(auto j = 0; j < col; j++){
c[i][j] = 0;
}
}
}
template <size_t row, size_t mid, size_t col>
void matmul1(float (&a)[row][mid], float (&b)[mid][col], float (&c)[row][col]){
for(auto i = 0; i < row; i++){
for(auto j = 0; j < col; j++){
float sum = 0;
for(auto k = 0; k < mid; k++){
sum += a[i][k] * b[k][j];
}
c[i][j] = sum;
}
}
}
template <size_t col>
inline void gebp(float (&Ab)[BLK_SIZE][BLK_SIZE], float (&Bp)[BLK_SIZE][col], float (&Cp)[BLK_SIZE][col]){
// We can optimize this subroutine but that'd be overkill for now.
for(auto j = 0; j < col; j++){
for(auto i = 0; i < BLK_SIZE; i++){
float sum = 0;
for(auto k = 0; k < BLK_SIZE; k++){
sum += Ab[i][k] * Bp[k][j];
}
Cp[i][j] += sum;
}
}
}
template <size_t row, size_t col>
inline void packb(float (&a)[row][col], float (&b)[BLK_SIZE][BLK_SIZE], size_t m, size_t n){
// size_t m, n in this case means the m,n-th block to pack.
auto start_row = m * BLK_SIZE;
auto start_col = n * BLK_SIZE;
for(size_t i = 0; i < BLK_SIZE; i++){
for(size_t j = 0; j < BLK_SIZE; j++){
b[i][j] = a[start_row + i][start_col + j];
}
}
}
template <size_t row, size_t mid, size_t col>
inline void gemm(float (&a)[row][mid], float (&b)[mid][col], float (&c)[row][col]){
// Divide up the matrix into panels:
// Suppose row / BLK_SIZE = M
// col / BLK_SIZE = N
// mid / BLK_SIZE = K
//
// For the rest of the function, we assume A, B, C as a, b, c variables,
// and {var}p = panel of var
// {var}b = block of var
//
// (TODO: We assume it's perfectly divisible for now)
// Layout: A = (M,K), B = (K,N), C = (M,N)
auto M = row / BLK_SIZE;
auto N = col / BLK_SIZE;
auto K = mid / BLK_SIZE;
for(auto p = 0; p < K; p++){
// Reassign B[p*BLK_SIZE : (p+1)*BLK_SIZE][:] into Bp
float (&Bp)[BLK_SIZE][col] = *(float (*)[BLK_SIZE][col]) &b[p * BLK_SIZE];
for(auto i = 0; i < M; i++){
// Pack A[i*BLK_SIZE : (i+1)*BLK_SIZE][p*BLK_SIZE : (p+1)*BLK_SIZE] into Ab
float Ab[BLK_SIZE][BLK_SIZE];
// Reassign C[i*BLK_SIZE : (i+1)*BLK_SIZE][:] into Cp
float (&Cp)[BLK_SIZE][col] = *(float (*)[BLK_SIZE][col]) &c[i * BLK_SIZE];
packb(a, Ab, i, p);
// The result of Ab and Bp should be in Cp
gebp(Ab, Bp, Cp);
}
}
}
template <size_t row, size_t col>
bool allclose(float (&a)[row][col], float (&b)[row][col], float threshold = 1e-5, bool verbose = true){
bool is_equal = true;
for(auto i = 0; i < row; i++){
for(auto j = 0; j < col; j++){
bool current_element = std::abs(a[i][j] - b[i][j]) < threshold;
if(verbose && !current_element){
std::cerr << "Element at [" << i << "][" << j << "] is incorrect : "
<< a[i][j] << " vs. " << b[i][j] << "." << std::endl;
}
is_equal = is_equal && current_element;
}
}
return is_equal;
}
float a[X_SIZE][Y_SIZE];
float b[Y_SIZE][Z_SIZE];
float c1[X_SIZE][Z_SIZE];
float c2[X_SIZE][Z_SIZE];
int main(){
initialize_matrices(a, b, c1);
TIME_IT(matmul1(a, b, c1))
// We must guarrantee c is all zeros at first.
initialize_matrices(a, b, c2);
TIME_IT(gemm(a, b, c2))
std::cout << allclose(c1, c2, 1e-1, true) << std::endl;
return 0;
}
g++-8 -o main main.cpp -std=c++17 -ffast-math -mavx2 -O3 -ftree-vectorize -fopt-info-vec-missed
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:1083:16: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/char_traits.h:320:25: note: not vectorized: not enough data-refs in basic block.
main.cpp:41:23: note: not vectorized: complicated access pattern.
main.cpp:41:23: note: bad data access.
main.cpp:42:27: note: Unknown misalignment, naturally aligned
main.cpp:36:23: note: not vectorized: complicated access pattern.
main.cpp:36:23: note: bad data access.
main.cpp:37:27: note: Unknown misalignment, naturally aligned
main.cpp:38:34: note: not vectorized: not enough data-refs in basic block.
main.cpp:38:13: note: not vectorized: no vectype for stmt: MEM[(float *)vectp_a.16_42] = vect__6.15_7;
scalar_type: vector(8) float
main.cpp:38:13: note: not vectorized: not enough data-refs in basic block.
main.cpp:36:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:43:34: note: not vectorized: not enough data-refs in basic block.
main.cpp:43:13: note: not vectorized: no vectype for stmt: MEM[(float *)vectp_b.10_15] = vect__12.9_38;
scalar_type: vector(8) float
main.cpp:43:13: note: not vectorized: not enough data-refs in basic block.
main.cpp:41:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:51:1: note: not vectorized: not enough data-refs in basic block.
main.cpp:121:23: note: not vectorized: multiple nested loops.
main.cpp:121:23: note: bad loop form.
main.cpp:125:27: note: not vectorized: multiple nested loops.
main.cpp:125:27: note: bad loop form.
main.cpp:69:23: note: not vectorized: multiple nested loops.
main.cpp:69:23: note: bad loop form.
main.cpp:70:27: note: step unknown.
main.cpp:70:27: note: not consecutive access *Cp_14[i_63][j_42] = _30;
main.cpp:70:27: note: not vectorized: complicated access pattern.
main.cpp:70:27: note: bad data access.
main.cpp:72:31: note: step unknown.
main.cpp:72:31: note: misalign = 0 bytes of ref Ab[i_63][k_64]
main.cpp:72:31: note: Unknown alignment for access: *Bp_12[k_64][j_42]
main.cpp:72:31: note: vector alignment may not be reachable
main.cpp:72:31: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:72:31: note: not ssa-name.
main.cpp:72:31: note: use not simple.
main.cpp:72:31: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:72:31: note: not ssa-name.
main.cpp:72:31: note: use not simple.
main.cpp:72:31: note: no array mode for V8SF[2048]
main.cpp:72:31: note: single-element interleaving not supported for not adjacent vector loads
main.cpp:72:31: note: not falling back to elementwise accesses
main.cpp:72:31: note: not vectorized: relevant stmt not supported: _24 = *Bp_12[k_64][j_42];
main.cpp:72:31: note: bad operation or unsupported loop bound.
main.cpp:72:31: note: step unknown.
main.cpp:72:31: note: misalign = 0 bytes of ref Ab[i_63][k_64]
main.cpp:72:31: note: Unknown alignment for access: *Bp_12[k_64][j_42]
main.cpp:72:31: note: vector alignment may not be reachable
main.cpp:72:31: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:72:31: note: not ssa-name.
main.cpp:72:31: note: use not simple.
main.cpp:72:31: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:72:31: note: not ssa-name.
main.cpp:72:31: note: use not simple.
main.cpp:72:31: note: no array mode for V4SF[2048]
main.cpp:72:31: note: single-element interleaving not supported for not adjacent vector loads
main.cpp:72:31: note: not falling back to elementwise accesses
main.cpp:72:31: note: not vectorized: relevant stmt not supported: _24 = *Bp_12[k_64][j_42];
main.cpp:72:31: note: bad operation or unsupported loop bound.
main.cpp:97:25: note: not vectorized: loop contains function calls or data references that cannot be analyzed
main.cpp:96:10: note: not vectorized: not enough data-refs in basic block.
main.cpp:95:10: note: not vectorized: not enough data-refs in basic block.
main.cpp:95:10: note: not vectorized: not enough data-refs in basic block.
main.cpp:97:25: note: not vectorized: not enough data-refs in basic block.
main.cpp:72:31: note: not consecutive access _24 = *Bp_12[k_64][j_42];
main.cpp:72:31: note: not consecutive access _23 = Ab[i_63][k_64];
main.cpp:72:31: note: not vectorized: no grouped stores in basic block.
main.cpp:70:27: note: not consecutive access _29 = *Cp_14[i_63][j_42];
main.cpp:70:27: note: not consecutive access *Cp_14[i_63][j_42] = _30;
main.cpp:70:27: note: not vectorized: no grouped stores in basic block.
main.cpp:69:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:125:27: note: not vectorized: not enough data-refs in basic block.
main.cpp:121:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:137:1: note: not vectorized: not enough data-refs in basic block.
main.cpp:140:6: note: not vectorized: control flow in loop.
main.cpp:140:6: note: bad loop form.
main.cpp:144:49: note: not vectorized: control flow in loop.
main.cpp:144:49: note: bad loop form.
main.cpp:145:13: note: not consecutive access _3 = *a_20(D)[i_17][j_131];
main.cpp:145:13: note: not consecutive access _4 = *b_21(D)[i_17][j_131];
main.cpp:145:13: note: not vectorized: no grouped stores in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:228:43: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:228:43: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/basic_ios.h:49:7: note: not consecutive access _118 = MEM[(const struct basic_ios *)_58]._M_ctype;
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/basic_ios.h:49:7: note: not consecutive access _55 = MEM[(struct basic_ostream *)_36]._vptr.basic_ostream;
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/basic_ios.h:49:7: note: not consecutive access _56 = MEM[(long int *)_55 + -24B];
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/basic_ios.h:49:7: note: not vectorized: no grouped stores in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:874:2: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:875:51: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:877:27: note: not consecutive access _127 = MEM[(const struct ctype *)_118].D.35451._vptr.facet;
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:877:27: note: not consecutive access _128 = MEM[(int (*) () *)_127 + 48B];
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:877:27: note: not vectorized: no grouped stores in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/bits/locale_facets.h:877:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:143:27: note: not vectorized: not enough data-refs in basic block.
main.cpp:142:23: note: not vectorized: not enough data-refs in basic block.
main.cpp:152:12: note: not vectorized: not enough data-refs in basic block.
main.cpp:55:23: note: not vectorized: multiple nested loops.
main.cpp:55:23: note: bad loop form.
main.cpp:56:27: note: step unknown.
main.cpp:56:27: note: inner step doesn't divide the vector alignment.
main.cpp:56:27: note: Unknown alignment for access: a[i_57][k_59]
main.cpp:56:27: note: misalign = 0 bytes of ref b[k_59][j_58]
main.cpp:56:27: note: misalign = 0 bytes of ref c1[i_57][j_58]
main.cpp:56:27: note: Unknown misalignment, naturally aligned
main.cpp:56:27: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:56:27: note: not ssa-name.
main.cpp:56:27: note: use not simple.
main.cpp:56:27: note: num. args = 4 (not unary/binary/ternary op).
main.cpp:56:27: note: not ssa-name.
main.cpp:56:27: note: use not simple.
main.cpp:56:27: note: can't use a fully-masked loop because no conditional operation is available.
main.cpp:162:5: note: not vectorized: not enough data-refs in basic block.
main.cpp:162:5: note: not vectorized: not enough data-refs in basic block.
main.cpp:58:31: note: not vectorized: no vectype for stmt: vect__37.134_72 = MEM[(float *)vectp_b.132_70];
scalar_type: vector(8) float
main.cpp:58:31: note: not consecutive access _36 = a[i_57][k_59];
main.cpp:58:31: note: not vectorized: no grouped stores in basic block.
main.cpp:61:13: note: not vectorized: no vectype for stmt: MEM[(float *)vectp_c1.138_80] = vect_sum_40.136_76;
scalar_type: vector(8) float
main.cpp:61:13: note: not vectorized: not enough data-refs in basic block.
main.cpp:55:23: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:562:44: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:561:18: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:562:44: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:561:18: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:175:29: note: not vectorized: not enough data-refs in basic block.
/usr/local/Cellar/gcc/8.2.0/include/c++/8.2.0/ostream:113:13: note: not vectorized: not enough data-refs in basic block.
main.cpp:170:1: note: not vectorized: not enough data-refs in basic block.
./main
Time elapsed : 4.03692
Time elapsed : 6.33257
我尝试使用gebp的块大小，看起来对于2048x1024 x 1024 x 2048矩阵，单个浮点似乎是-O2和没有矢量化指令的最佳性能。这有点奇怪，因为我认为更大的块可以驻留在tlb / L1 / L2缓存中。
1 float:
Time elapsed : 2.33989
8 floats:
Time elapsed : 4.34408
16 floats:
Time elapsed : 4.522
32 floats:
Time elapsed : 6.42583
64 floats:
Time elapsed : 7.35193
128 floats:
Time elapsed : 8.18502
256 floats:
Time elapsed : 8.1686
但是，这里有一些非常有趣的东西，带有矢量化指令：
1 float:
Time elapsed : 1.3284
2 floats:
Time elapsed : 1.62332
4 floats:
Time elapsed : 0.421444
8 floats:
Time elapsed : 0.527425
16 floats:
Time elapsed : 4.03877
......从这里变得更糟似乎gcc自动矢量化器似乎已针对AVX和AVX2指令的特殊情况4和8进行了积极优化。这只是一个推测，而不是一个严谨的发现。在极端情况下，矢量化指令的加速实际上从4倍到8倍！
在矢量化环境下的2048x2048方阵乘法：
Time elapsed : 11.9983 <-- naive
Time elapsed : 0.822692 <-- vectorized
Time elapsed : 0.472855 <-- ATLAS blas GEMM
这是一个惊人的> 10倍的速度提升。与ATLAS相比，在我的基准测试中（使用xtensor作为接口），我的实现速度仅慢了2倍。
只是一个在评论中看起来不太好的假设：
auto start_row = m * BLK_SIZE;
auto start_col = n * BLK_SIZE;
for(size_t i = 0; i < BLK_SIZE; i++){
for(size_t j = 0; j < BLK_SIZE; j++){
b[i][j] = a[start_row + i][start_col + j];
}
}
可以略微优化成
const auto start_row = m * BLK_SIZE;
const auto start_col = n * BLK_SIZE;
const auto end_row = start_row + BLK_SIZE;
const auto end_col = start_col + BLK_SIZE;
for(size_t i = start_row; i < end_row; i++){
for(size_t j = start_col; j < end_col; j++){
b[i][j] = a[i][j];
}
}
它可能已经由编译器完成但谁知道没有检查汇编输出？
还不是memcpy好在这里？它看起来像一个副本。