为什么幼稚的C ++矩阵乘法比BLAS慢100倍？
我正在研究大矩阵乘法，并进行以下实验以形成基线测试：
以下是天真的C ++实现：
#include <iostream>
#include <algorithm>
using namespace std;
int main()
{
constexpr size_t dim = 4096;
float* x = new float[dim*dim];
float* y = new float[dim*dim];
float* z = new float[dim*dim];
random_device rd;
mt19937 gen(rd());
normal_distribution<float> dist(0, 1);
for (size_t i = 0; i < dim*dim; i++)
{
x[i] = dist(gen);
y[i] = dist(gen);
}
for (size_t row = 0; row < dim; row++)
for (size_t col = 0; col < dim; col++)
{
float acc = 0;
for (size_t k = 0; k < dim; k++)
acc += x[row*dim + k] * y[k*dim + col];
z[row*dim + col] = acc;
}
float t = 0;
for (size_t i = 0; i < dim*dim; i++)
t += z[i];
cout << t << endl;
delete x;
delete y;
delete z;
}
编译并运行：
$ g++ -std=gnu++11 -O3 test.cpp
$ time ./a.out
这是Octave / matlab实现：
X = stdnormal_rnd(4096, 4096);
Y = stdnormal_rnd(4096, 4096);
Z = X*Y;
sum(sum(Z))
跑：
$ time octave < test.octave
引擎盖下的八度使用BLAS（我假设sgemm功能）
硬件是Linux x86-64上配备24 GB RAM的i7 3930X。 BLAS似乎在使用两个内核。 也许是超线程对？
我发现使用-O3上的GCC 4.7编译的C ++版本需要9分钟才能执行：
real 9m2.126s
user 9m0.302s
sys 0m0.052s
八度版本花了6秒钟：
real 0m5.985s
user 0m10.881s
sys 0m0.144s
据我所知，BLAS已经对所有地狱进行了优化，而天真的算法完全忽略了缓存等等，但是认真 - 90次？
任何人都可以解释这种差异？ BLAS实现的结构究竟是什么？ 我看到它正在使用Fortran，但在CPU级别发生了什么？ 它使用什么算法？ 它如何使用CPU缓存？ 它调用什么x86-64机器指令？ （它是否使用像AVX这样的高级CPU功能？）它从哪里获得这个额外的速度？
对C ++算法的哪些关键优化可以使其与BLAS版本保持一致？
我在gdb下运行了八度，并通过几次计算停止了它。 它已经开始了第二个线程，这里是堆栈（所有的停顿看起来都很相似）：
(gdb) thread 1
#0 0x00007ffff6e17148 in pthread_join () from /lib/x86_64-linux-gnu/libpthread.so.0
#1 0x00007ffff1626721 in ATL_join_tree () from /usr/lib/libblas.so.3
#2 0x00007ffff1626702 in ATL_join_tree () from /usr/lib/libblas.so.3
#3 0x00007ffff15ae357 in ATL_dptgemm () from /usr/lib/libblas.so.3
#4 0x00007ffff1384b59 in atl_f77wrap_dgemm_ () from /usr/lib/libblas.so.3
#5 0x00007ffff193effa in dgemm_ () from /usr/lib/libblas.so.3
#6 0x00007ffff6049727 in xgemm(Matrix const&, Matrix const&, blas_trans_type, blas_trans_type) () from /usr/lib/x86_64-linux-gnu/liboctave.so.1
#7 0x00007ffff6049954 in operator*(Matrix const&, Matrix const&) () from /usr/lib/x86_64-linux-gnu/liboctave.so.1
#8 0x00007ffff7839e4e in ?? () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#9 0x00007ffff765a93a in do_binary_op(octave_value::binary_op, octave_value const&, octave_value const&) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#10 0x00007ffff76c4190 in tree_binary_expression::rvalue1(int) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#11 0x00007ffff76c33a5 in tree_simple_assignment::rvalue1(int) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#12 0x00007ffff76d0864 in tree_evaluator::visit_statement(tree_statement&) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#13 0x00007ffff76cffae in tree_evaluator::visit_statement_list(tree_statement_list&) () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#14 0x00007ffff757f6d4 in main_loop() () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
#15 0x00007ffff7527abf in octave_main () from /usr/lib/x86_64-linux-gnu/liboctinterp.so.1
(gdb) thread 2
#0 0x00007ffff14ba4df in ATL_dJIK56x56x56TN56x56x0_a1_b1 () from /usr/lib/libblas.so.3
(gdb) bt
#0 0x00007ffff14ba4df in ATL_dJIK56x56x56TN56x56x0_a1_b1 () from /usr/lib/libblas.so.3
#1 0x00007ffff15a5fd7 in ATL_dmmIJK2 () from /usr/lib/libblas.so.3
#2 0x00007ffff15a6ae4 in ATL_dmmIJK () from /usr/lib/libblas.so.3
#3 0x00007ffff1518e65 in ATL_dgemm () from /usr/lib/libblas.so.3
#4 0x00007ffff15adf7a in ATL_dptgemm0 () from /usr/lib/libblas.so.3
#5 0x00007ffff6e15e9a in start_thread () from /lib/x86_64-linux-gnu/libpthread.so.0
#6 0x00007ffff6b41cbd in clone () from /lib/x86_64-linux-gnu/libc.so.6
#7 0x0000000000000000 in ?? ()
正如预期的那样，它正在呼唤BLAS gemm 。
第一个线程似乎加入了第二个线程，所以我不确定这两个线程是否占用了200％的CPU使用率。
哪个库是ATL_dgemm libblas.so.3，它的代码在哪里？
$ ls -al /usr/lib/libblas.so.3
/usr/lib/libblas.so.3 -> /etc/alternatives/libblas.so.3
$ ls -al /etc/alternatives/libblas.so.3
/etc/alternatives/libblas.so.3 -> /usr/lib/atlas-base/atlas/libblas.so.3
$ ls -al /usr/lib/atlas-base/atlas/libblas.so.3
/usr/lib/atlas-base/atlas/libblas.so.3 -> libblas.so.3.0
$ ls -al /usr/lib/atlas-base/atlas/libblas.so.3.0
/usr/lib/atlas-base/atlas/libblas.so.3.0
$ dpkg -S /usr/lib/atlas-base/atlas/libblas.so.3.0
libatlas3-base: /usr/lib/atlas-base/atlas/libblas.so.3.0
$ apt-get source libatlas3-base
这是ATLAS 3.8.4
以下是我后来实现的优化：
使用平铺方法将X，Y和Z的64x64块预加载到单独的阵列中。
更改每个块的计算，以便内部循环如下所示：
for (size_t tcol = 0; tcol < block_width; tcol++)
bufz[trow][tcol] += B * bufy[tk][tcol];
这允许GCC自动向量化SIMD指令，并允许指令级并行（我认为）。
开启march=corei7-avx 。 这增加了30％的额外速度，但是作弊，因为我认为BLAS库是预先构建的。
代码如下：
#include <iostream>
#include <algorithm>
using namespace std;
constexpr size_t dim = 4096;
constexpr size_t block_width = 64;
constexpr size_t num_blocks = dim / block_width;
double X[dim][dim], Y[dim][dim], Z[dim][dim];
double bufx[block_width][block_width];
double bufy[block_width][block_width];
double bufz[block_width][block_width];
void calc_block()
{
for (size_t trow = 0; trow < block_width; trow++)
for (size_t tk = 0; tk < block_width; tk++)
{
double B = bufx[trow][tk];
for (size_t tcol = 0; tcol < block_width; tcol++)
bufz[trow][tcol] += B * bufy[tk][tcol];
}
}
int main()
{
random_device rd;
mt19937 gen(rd());
normal_distribution<double> dist(0, 1);
for (size_t row = 0; row < dim; row++)
for (size_t col = 0; col < dim; col++)
{
X[row][col] = dist(gen);
Y[row][col] = dist(gen);
Z[row][col] = 0;
}
for (size_t block_row = 0; block_row < num_blocks; block_row++)
for (size_t block_col = 0; block_col < num_blocks; block_col++)
{
for (size_t trow = 0; trow < block_width; trow++)
for (size_t tcol = 0; tcol < block_width; tcol++)
bufz[trow][tcol] = 0;
for (size_t block_k = 0; block_k < num_blocks; block_k++)
{
for (size_t trow = 0; trow < block_width; trow++)
for (size_t tcol = 0; tcol < block_width; tcol++)
{
bufx[trow][tcol] = X[block_row*block_width + trow][block_k*block_width + tcol];
bufy[trow][tcol] = Y[block_k*block_width + trow][block_col*block_width + tcol];
}
calc_block();
}
for (size_t trow = 0; trow < block_width; trow++)
for (size_t tcol = 0; tcol < block_width; tcol++)
Z[block_row*block_width + trow][block_col*block_width + tcol] = bufz[trow][tcol];
}
double t = 0;
for (size_t row = 0; row < dim; row++)
for (size_t col = 0; col < dim; col++)
t += Z[row][col];
cout << t << endl;
}
所有的动作都在calc_block函数中 - 超过90％的时间花在了它上面。
新的时间是：
real 0m17.370s
user 0m17.213s
sys 0m0.092s
这更接近。
calc_block函数的反编译如下：
0000000000401460 <_Z10calc_blockv>:
401460: b8 e0 21 60 00 mov $0x6021e0,%eax
401465: 41 b8 e0 23 61 00 mov $0x6123e0,%r8d
40146b: 31 ff xor %edi,%edi
40146d: 49 29 c0 sub %rax,%r8
401470: 49 8d 34 00 lea (%r8,%rax,1),%rsi
401474: 48 89 f9 mov %rdi,%rcx
401477: ba e0 a1 60 00 mov $0x60a1e0,%edx
40147c: 48 c1 e1 09 shl $0x9,%rcx
401480: 48 81 c1 e0 21 61 00 add $0x6121e0,%rcx
401487: 66 0f 1f 84 00 00 00 nopw 0x0(%rax,%rax,1)
40148e: 00 00
401490: c4 e2 7d 19 01 vbroadcastsd (%rcx),%ymm0
401495: 48 83 c1 08 add $0x8,%rcx
401499: c5 fd 59 0a vmulpd (%rdx),%ymm0,%ymm1
40149d: c5 f5 58 08 vaddpd (%rax),%ymm1,%ymm1
4014a1: c5 fd 29 08 vmovapd %ymm1,(%rax)
4014a5: c5 fd 59 4a 20 vmulpd 0x20(%rdx),%ymm0,%ymm1
4014aa: c5 f5 58 48 20 vaddpd 0x20(%rax),%ymm1,%ymm1
4014af: c5 fd 29 48 20 vmovapd %ymm1,0x20(%rax)
4014b4: c5 fd 59 4a 40 vmulpd 0x40(%rdx),%ymm0,%ymm1
4014b9: c5 f5 58 48 40 vaddpd 0x40(%rax),%ymm1,%ymm1
4014be: c5 fd 29 48 40 vmovapd %ymm1,0x40(%rax)
4014c3: c5 fd 59 4a 60 vmulpd 0x60(%rdx),%ymm0,%ymm1
4014c8: c5 f5 58 48 60 vaddpd 0x60(%rax),%ymm1,%ymm1
4014cd: c5 fd 29 48 60 vmovapd %ymm1,0x60(%rax)
4014d2: c5 fd 59 8a 80 00 00 vmulpd 0x80(%rdx),%ymm0,%ymm1
4014d9: 00
4014da: c5 f5 58 88 80 00 00 vaddpd 0x80(%rax),%ymm1,%ymm1
4014e1: 00
4014e2: c5 fd 29 88 80 00 00 vmovapd %ymm1,0x80(%rax)
4014e9: 00
4014ea: c5 fd 59 8a a0 00 00 vmulpd 0xa0(%rdx),%ymm0,%ymm1
4014f1: 00
4014f2: c5 f5 58 88 a0 00 00 vaddpd 0xa0(%rax),%ymm1,%ymm1
4014f9: 00
4014fa: c5 fd 29 88 a0 00 00 vmovapd %ymm1,0xa0(%rax)
401501: 00
401502: c5 fd 59 8a c0 00 00 vmulpd 0xc0(%rdx),%ymm0,%ymm1
401509: 00
40150a: c5 f5 58 88 c0 00 00 vaddpd 0xc0(%rax),%ymm1,%ymm1
401511: 00
401512: c5 fd 29 88 c0 00 00 vmovapd %ymm1,0xc0(%rax)
401519: 00
40151a: c5 fd 59 8a e0 00 00 vmulpd 0xe0(%rdx),%ymm0,%ymm1
401521: 00
401522: c5 f5 58 88 e0 00 00 vaddpd 0xe0(%rax),%ymm1,%ymm1
401529: 00
40152a: c5 fd 29 88 e0 00 00 vmovapd %ymm1,0xe0(%rax)
401531: 00
401532: c5 fd 59 8a 00 01 00 vmulpd 0x100(%rdx),%ymm0,%ymm1
401539: 00
40153a: c5 f5 58 88 00 01 00 vaddpd 0x100(%rax),%ymm1,%ymm1
401541: 00
401542: c5 fd 29 88 00 01 00 vmovapd %ymm1,0x100(%rax)
401549: 00
40154a: c5 fd 59 8a 20 01 00 vmulpd 0x120(%rdx),%ymm0,%ymm1
401551: 00
401552: c5 f5 58 88 20 01 00 vaddpd 0x120(%rax),%ymm1,%ymm1
401559: 00
40155a: c5 fd 29 88 20 01 00 vmovapd %ymm1,0x120(%rax)
401561: 00
401562: c5 fd 59 8a 40 01 00 vmulpd 0x140(%rdx),%ymm0,%ymm1
401569: 00
40156a: c5 f5 58 88 40 01 00 vaddpd 0x140(%rax),%ymm1,%ymm1
401571: 00
401572: c5 fd 29 88 40 01 00 vmovapd %ymm1,0x140(%rax)
401579: 00
40157a: c5 fd 59 8a 60 01 00 vmulpd 0x160(%rdx),%ymm0,%ymm1
401581: 00
401582: c5 f5 58 88 60 01 00 vaddpd 0x160(%rax),%ymm1,%ymm1
401589: 00
40158a: c5 fd 29 88 60 01 00 vmovapd %ymm1,0x160(%rax)
401591: 00
401592: c5 fd 59 8a 80 01 00 vmulpd 0x180(%rdx),%ymm0,%ymm1
401599: 00
40159a: c5 f5 58 88 80 01 00 vaddpd 0x180(%rax),%ymm1,%ymm1
4015a1: 00
4015a2: c5 fd 29 88 80 01 00 vmovapd %ymm1,0x180(%rax)
4015a9: 00
4015aa: c5 fd 59 8a a0 01 00 vmulpd 0x1a0(%rdx),%ymm0,%ymm1
4015b1: 00
4015b2: c5 f5 58 88 a0 01 00 vaddpd 0x1a0(%rax),%ymm1,%ymm1
4015b9: 00
4015ba: c5 fd 29 88 a0 01 00 vmovapd %ymm1,0x1a0(%rax)
4015c1: 00
4015c2: c5 fd 59 8a c0 01 00 vmulpd 0x1c0(%rdx),%ymm0,%ymm1
4015c9: 00
4015ca: c5 f5 58 88 c0 01 00 vaddpd 0x1c0(%rax),%ymm1,%ymm1
4015d1: 00
4015d2: c5 fd 29 88 c0 01 00 vmovapd %ymm1,0x1c0(%rax)
4015d9: 00
4015da: c5 fd 59 82 e0 01 00 vmulpd 0x1e0(%rdx),%ymm0,%ymm0
4015e1: 00
4015e2: c5 fd 58 80 e0 01 00 vaddpd 0x1e0(%rax),%ymm0,%ymm0
4015e9: 00
4015ea: 48 81 c2 00 02 00 00 add $0x200,%rdx
4015f1: 48 39 ce cmp %rcx,%rsi
4015f4: c5 fd 29 80 e0 01 00 vmovapd %ymm0,0x1e0(%rax)
4015fb: 00
4015fc: 0f 85 8e fe ff ff jne 401490 <_Z10calc_blockv+0x30>
401602: 48 83 c7 01 add $0x1,%rdi
401606: 48 05 00 02 00 00 add $0x200,%rax
40160c: 48 83 ff 40 cmp $0x40,%rdi
401610: 0f 85 5a fe ff ff jne 401470 <_Z10calc_blockv+0x10>
401616: c5 f8 77 vzeroupper
401619: c3 retq
40161a: 66 0f 1f 44 00 00 nopw 0x0(%rax,%rax,1)
以下是导致代码和BLAS之间性能差异的三个因素（加上Strassen算法的注释）。
在你的内循环中，在k ，你有y[k*dim + col] 。 由于内存缓存的排列方式， k连续值与相同的缓存集具有相同的dim和col映射。 缓存结构的方式是，每个内存地址都有一个缓存集，缓存中的内容必须保存。 每个缓存集有几行（四个是典型的数字），并且每行都可以保存任何映射到该特定缓存集的内存地址。
因为内部循环以这种方式迭代y ，所以每次它使用y一个元素时，它必须将该元素的内存加载到与上一次迭代相同的集合中。 这迫使该集合中之前的一个缓存行被驱逐。 然后，在col循环的下一次迭代中， y所有元素都已从缓存中逐出，因此它们必须重新加载。
因此， 每次你的循环加载一个y元素时，它都必须从内存中加载，这需要很多CPU周期。
高性能代码通过两种方式避免了这一点。 其一，它将工作分成更小的块。 行和列被分割成更小的大小，并用较短的循环进行处理，这些循环能够使用高速缓存行中的所有元素，并在每个元素进入下一个块之前多次使用每个元素。 因此，对x元素和y元素的大多数引用来自缓存，通常在单个处理器周期中。 二，在某些情况下，代码会将数据从矩阵的一列（由于几何形状而导致高速缓存抖动）复制到临时缓冲区的一行（这可避免抖动）。 这再次允许大多数内存引用从缓存而不是从内存提供。
另一个因素是使用单指令多数据（SIMD）功能。 许多现代处理器都有指令，可以在一条指令中加载多个元素（四个float元素是典型的，但现在有八个），存储多个元素，添加多个元素（例如，对于这四个元素中的每一个，将其添加到相应的一个元素四个），乘以多个元素，等等。 只要使用这些说明，即可让您的代码快速四倍，前提是您可以安排使用这些说明。
这些指令不能在标准C中直接访问。现在，一些优化器尽可能地使用这些指令，但这种优化很困难，并且从它中获益不大。 许多编译器提供了允许访问这些指令的语言的扩展。 就个人而言，我通常更喜欢用汇编语言编写以使用SIMD。
另一个因素是在处理器上使用指令级并行执行功能。 注意到在你的内部循环中， acc被更新了。 下一次迭代不能添加到acc ，直到前一次迭代已经完成更新acc 。 高性能代码将保持并行运行多个和（甚至多个SIMD和）。 这样做的结果是，当执行一个总和的加法时，将开始另一个总和的加法。 在当今的处理器中，一次支持四个或更多的浮点操作是很常见的。 正如所写，你的代码根本无法做到这一点。 一些编译器会尝试通过重新排列循环来优化代码，但是这需要编译器能够看到特定循环的迭代彼此独立，或者可以通过另一个循环等等来减少。
使用高速缓存有效地提供了十倍的性能改进是非常可行的，SIMD提供了另外四个，指令级并行性提供了另外四个，共提供了160个。
这是基于维基百科页面对Strassen算法效果的一个非常粗略的估计。 维基百科页面说Strassen比n = 100时的直接乘法略好。这表明执行时间常数因子的比率是1003 /1002.807≈2.4。 很明显，这取决于处理器型号，与缓存效果交互的矩阵大小等等，会有很大的差异。 然而，简单的推断表明Strassen在n = 4096（（4096/100）3-2.807≈2.05）时的直接乘法性能是其两倍。 再次，这只是一个球估计。
至于后面的优化，请考虑内部循环中的这段代码：
bufz[trow][tcol] += B * bufy[tk][tcol];
与此相关的一个潜在问题是， bufz可能重叠bufy 。 既然你使用了bufz和bufy全局定义，编译器可能知道它们在这种情况下不重叠。 但是，如果将此代码移到传递bufz和bufy作为参数的子例程中，并且特别是在单独的源文件中编译该子例程时，编译器不太可能知道bufz和bufy不重叠。 在这种情况下，编译器无法矢量化或重新排序代码，因为此迭代中的bufz[trow][tcol]可能与另一次迭代中的bufy[tk][tcol]相同。
即使编译器可以看到子程序在当前源模块中被调用了不同的bufz和bufy ，但如果该子程序具有extern连接（缺省），则编译器必须允许从外部模块调用子程序，因此如果bufz和bufy重叠，它必须生成正确工作的代码。 （编译器可以处理的一种方式是生成两个版本的例程，一个从外部模块调用，另一个从当前正在编译的模块调用。无论是否依赖于编译器，优化开关，等等）。如果你声明这个例程是static ，那么编译器知道它不能从外部模块调用（除非你使用它的地址，并且有可能地址在当前模块之外传递）。
另一个潜在的问题是，即使编译器将此代码向量化，也不一定会为您执行的处理器生成最佳代码。 查看生成的汇编代码，看起来编译器只反复使用%ymm1 。 一遍又一遍，它把从内存的值到%ymm1 ，增加了从内存的值%ymm1 ，并存储从价值%ymm1内存。 这有两个问题。
一，你不希望这些部分总和经常存储到内存中。 你需要在寄存器中累加很多新的内容，而寄存器很少会被写入内存。 说服编译器做到这一点可能需要重写代码，以明确在循环完成后将部分和保存在临时对象中并将它们写入内存。
二，这些指令名义上是连续相关的。 直到乘法完成后才能开始添加，直到添加完成后，存储才能写入内存。 Core i7具有很强的无序执行能力。 所以，当它增加了等待开始执行时，它会在指令流中稍后查看乘法并启动它。 （即使这个乘法也使用%ymm1 ，处理器会立即重新映射寄存器，以便它使用不同的内部寄存器来执行此乘法操作。）即使您的代码被连续的依赖项填充，处理器也会尝试执行几条指令立刻。 然而，一些事情可能会干扰这一点。 您可以用尽处理器用于重命名的内部寄存器。 您使用的内存地址可能会遇到错误冲突。 （处理器查看十几个内存地址的低位，以查看该地址是否与另一个尝试从较早的指令加载或存储的地址相同。如果这些位相同，则处理器具有以延迟当前的加载或存储，直到它可以验证整个地址是不同的，这种延迟可能会超过当前的加载或存储）。因此，最好是有明显独立的指令。
这是我更愿意在汇编中编写高性能代码的另一个原因。 要用C语言来完成，你必须说服编译器给你这样的指示，通过编写你自己的一些SIMD代码（使用它们的语言扩展）和手动展开循环（写出多次迭代）。
复制进入和退出缓冲区时，可能会出现类似的问题。 但是，您报告90％的时间都花在calc_block ，所以我没有仔细看过。
另外，斯特拉森的算法是否解释了其余的差异？
Strassen的算法比naïve算法有两个优点：
B*M½ ，其中B是高速缓存行大小，M是高速缓存大小。 我认为第二点对于你正在经历的放缓有很大的影响。 如果你在Linux下运行你的应用程序，我建议你使用perf工具运行它们，它告诉你程序遇到了多少缓存未命中。
我不知道这些信息有多可靠，但维基百科说BLAS使用Strassen算法来处理大矩阵。 你的确很大。 那大约是O（n ^ 2.807），比你的O（n ^ 3）天真的算法要好。
链接地址: http://www.djcxy.com/p/15081.html上一篇: Why is a naïve C++ matrix multiplication 100 times slower than BLAS?
下一篇: Why does the speed of memcpy() drop dramatically every 4KB?