矩阵乘法的速度取决于愚蠢的事情
matrix multiplication speed depends on silly things
我写了一个快速矩阵乘法的程序。为了最大限度地使用CPU缓存，它将矩阵划分为10*10个单元，并分别乘以每个单元(将单元大小增加到20*20或50*50不会显著改变时间)。
结果表明，速度在很大程度上取决于是否提前知道矩阵大小和单元大小。
程序是：
#include <cmath>
#include <cstdlib>
#include <iostream>
using namespace std;
#define forall(i,n) for(int i=0; i<(int)(n); i++)
inline void Load(int N, int N2, float* x2, float* x, int iStart, int jStart) {
int start = iStart * N + jStart;
forall (i, N2)
forall (j, N2)
x2[i * N2 + j] = x[start + i * N + j];
}
inline void Add(int N, int N2, float* x, float* x2, int iStart, int jStart) {
int start = iStart * N + jStart;
forall (i, N2)
forall (j, N2)
x[start + i * N + j] += x2[i * N2 + j];
}
inline void Mul(int N, float* z, float* x, float* y) {
forall (i, N)
forall (j, N) {
double sum = 0;
forall (k, N)
sum += x[i*N+k] * y[k*N+j];
z[i*N+j] = sum;
}
}
inline double RandReal() {return random()/((double)RAND_MAX+1);}
int main(int argc, char** argv) {
#if VAR==3
const int N = atoi(argv[1]), N2 = atoi(argv[2]);
#elif VAR==2
const int N = 2000, N2 = atoi(argv[2]);
#elif VAR==1
const int N = atoi(argv[1]), N2 = 10;
#elif VAR==0
const int N = 2000, N2 = 10;
#else
#error "Bad VAR"
#endif
cout << "VAR=" << VAR << " N=" << N << " N2=" << N2 << endl;
float x[N*N], y[N*N];
forall (i, N)
forall (j, N) {
x[i*N+j] = RandReal();
y[i*N+j] = RandReal();
}
float z[N*N];
forall (i, N)
forall (j, N)
z[i*N+j] = 0;
for (int i1 = 0; i1 < N; i1 += N2) {
float x2[N2*N2], y2[N2*N2], z2[N2*N2];
for (int j1 = 0; j1 < N; j1 += N2) {
Load(N, N2, x2, x, i1, j1);
for (int k1 = 0; k1 < N; k1 += N2) {
Load(N, N2, y2, y, j1, k1);
Mul(N2, z2, x2, y2);
Add(N, N2, z, z2, i1, k1);
}
}
}
double val = 0, val2 = 0;
forall (i, N)
forall (j, N)
val += z[i*N+j], val2 += z[i*N+j]*(i+j);
cout << "val=" << val << " val2=" << val2 << endl;
}
现在执行时间：
$ for a in 0 1 2 3; do g++ -DVAR=$a -O3 -Wall -o mat mat.cpp; time ./mat 2000 10; done
VAR=0 N=2000 N2=10
val=2.00039e+09 val2=3.99867e+12
real 0m8.127s
user 0m8.108s
sys 0m0.020s
VAR=1 N=2000 N2=10
val=2.00039e+09 val2=3.99867e+12
real 0m3.304s
user 0m3.292s
sys 0m0.012s
VAR=2 N=2000 N2=10
val=2.00039e+09 val2=3.99867e+12
real 0m25.395s
user 0m25.388s
sys 0m0.008s
VAR=3 N=2000 N2=10
val=2.00039e+09 val2=3.99867e+12
real 0m25.515s
user 0m25.495s
sys 0m0.016s
简单来说：
- 不知道矩阵大小，知道单元格大小：3.3秒
- 了解矩阵大小和单元格大小：8.1秒
- 不知道单元格大小：25.5秒
为什么？我使用g++5.4.0.
inline不起任何作用，如果我们去掉它，结果会是一样的。
引言：这篇文章的大部分内容都被重写了，所以下面的一些评论不再有多大意义了。如果你介意的话，请查看编辑后的详细信息。所以…
tl；dr
- 配置文件以查找热循环
- 如果可能的话，使用常量计数
- 如果没有，请尝试手动展开它们
- OP的结果是个谜，没有人能复制出如此极端的东西
我同意@user4581301的观点——编译器提前知道的越多，就越能为您优化代码。
但你需要了解一下这个代码——墙上的时钟时间只会带你走这么远。我对gcc的profiler了解不多(我有一个很好的用于MSVC的profiler)，但你可以在这里碰碰运气。
使用Godbolt作为一种工具，尝试学习一些汇编程序也是值得的(正如@RetiredNinja一开始所说)，尤其是如果你想理解像这样戏剧性的放缓。
说了这么多，你的时间安排对我来说毫无意义，所以你的系统正在发生一些奇怪的事情。所以我自己做了一些测试，结果和你的明显不同。我在MSVC上运行了其中一些测试(因为我在那里有非常棒的评测工具)，在Mac上的gcc上运行了一些测试(尽管我认为这实际上是秘密的)。对不起，我没有linux盒子。
首先，让我们来处理在堆栈上分配如此大的对象的问题。这显然是不明智的，我在MSVC上根本无法做到这一点，因为它不支持可变长度数组，但我在Mac上的测试表明，一旦我增加了ulimit以使其工作，这对时间没有任何影响(请参阅此处)。所以我认为我们可以忽略这一因素，正如你自己在评论中所说的那样。
现在让我们看看我在Mac上的时间：
VAR=0 USE_STACK=0 N=2000 (known) N2=10 (known)
user 0m10.813s
VAR=1 USE_STACK=0 N=2000 (unknown) N2=10 (known)
user 0m11.008s
VAR=2 USE_STACK=0 N=2000 (known) N2=10 (unknown)
user 0m12.714s
VAR=3 USE_STACK=0 N=2000 (unknown) N2=10 (unknown)
user 0m12.778s
VAR=0 USE_STACK=1 N=2000 (known) N2=10 (known)
user 0m10.617s
VAR=1 USE_STACK=1 N=2000 (unknown) N2=10 (known)
user 0m10.987s
VAR=2 USE_STACK=1 N=2000 (known) N2=10 (unknown)
user 0m12.653s
VAR=3 USE_STACK=1 N=2000 (unknown) N2=10 (unknown)
user 0m12.673s
那里没什么可看的；让我们继续看我在MSVC上观察到的内容(在那里我可以评测)：
VAR=0 USE_STACK=0 N=2000 (known) N2=10 (known)
Elapsed: 0:00:06.89
VAR=1 USE_STACK=0 N=2000 (unknown) N2=10 (known)
Elapsed: 0:00:06.86
VAR=2 USE_STACK=0 N=2000 (known) N2=10 (unknown)
Elapsed: 0:00:10.24
VAR=3 USE_STACK=0 N=2000 (unknown) N2=10 (unknown)
Elapsed: 0:00:10.39
现在我们有一些东西可以咬进去。正如@geza所观察到的，当N2未知时，代码需要更长的时间才能运行，这完全符合人们的预期，因为这是热循环的位置，而且当编译器知道它实际上由少量已知迭代组成时，它更有可能展开这样的循环。
因此，让我们从探查器中获取一些信息。它告诉我，热循环(相当长的一段时间)是Mul():中的内循环
inline void Mul(int N, float* z, float* x, float* y) {
forall (i, N)
forall (j, N) {
double sum = 0;
=>对于所有(k，N)=>sum+=x[i*N+k]*y[kN+j]；z[iN+j]=总和；}}
同样，我不能说这让我很惊讶，当我查看代码时，我可以看到循环根本没有展开(为了简洁起见，省略了循环设置代码)：
$1:
movss xmm0,dword ptr [r9+rsi*4]
mulss xmm0,dword ptr [r8+4]
movss xmm1,dword ptr [r9+r15*4]
mulss xmm1,dword ptr [r8]
cvtps2pd xmm2,xmm0
cvtps2pd xmm0,xmm1
movss xmm1,dword ptr [r8+8]
mulss xmm1,dword ptr [r9]
addsd xmm0,xmm3
addsd xmm2,xmm0
cvtps2pd xmm0,xmm1
movss xmm1,dword ptr [r9+r14*4]
movaps xmm3,xmm2
mulss xmm1,dword ptr [r8+0Ch]
add r9,rbp
add r8,10h
addsd xmm3,xmm0
cvtps2pd xmm0,xmm1
addsd xmm3,xmm0
sub rcx,1
jne $1
现在看来，展开该循环不会节省任何费用，因为与执行其中的所有其他代码相比，循环将是廉价的，但如果你看看N2已知时对同一循环的反汇编，你会发现一个惊喜：
movss xmm0,dword ptr [rax-8]
mulss xmm0,dword ptr [rcx-50h]
cvtps2pd xmm2,xmm0
movss xmm0,dword ptr [rcx-28h]
mulss xmm0,dword ptr [rax-4]
addsd xmm2,xmm7
cvtps2pd xmm1,xmm0
movss xmm0,dword ptr [rcx]
mulss xmm0,dword ptr [rax]
addsd xmm2,xmm1
cvtps2pd xmm1,xmm0
movss xmm0,dword ptr [rcx+28h]
mulss xmm0,dword ptr [rax+4]
addsd xmm2,xmm1
cvtps2pd xmm1,xmm0
movss xmm0,dword ptr [rcx+50h]
mulss xmm0,dword ptr [rax+8]
addsd xmm2,xmm1
cvtps2pd xmm1,xmm0
movss xmm0,dword ptr [rcx+78h]
mulss xmm0,dword ptr [rax+0Ch]
addsd xmm2,xmm1
cvtps2pd xmm1,xmm0
movss xmm0,dword ptr [rcx+0A0h]
mulss xmm0,dword ptr [rax+10h]
addsd xmm2,xmm1
cvtps2pd xmm1,xmm0
movss xmm0,dword ptr [rcx+0C8h]
mulss xmm0,dword ptr [rax+14h]
addsd xmm2,xmm1
cvtps2pd xmm1,xmm0
movss xmm0,dword ptr [rcx+0F0h]
mulss xmm0,dword ptr [rax+18h]
addsd xmm2,xmm1
cvtps2pd xmm1,xmm0
movss xmm0,dword ptr [rcx+118h]
mulss xmm0,dword ptr [rax+1Ch]
addsd xmm2,xmm1
cvtps2pd xmm1,xmm0
addsd xmm2,xmm1
cvtpd2ps xmm0,xmm2
movss dword ptr [rdx+rcx],xmm0
现在没有循环，将要执行的指令数量明显减少。也许MS毕竟不是一群愚蠢的傻瓜。
最后，作为练习，让我们手动快速展开循环，看看我们得到了什么时间(我没有查看生成的代码)：
#define UNROLL_LOOP 1
inline void Mul(int N, float* z, float* x, float* y) {
forall (i, N)
forall (j, N) {
double sum = 0;
#if UNROLL_LOOP
assert (N == 10);
sum += x[i*N] * y[0*N+j];
sum += x[i*N+1] * y[1*N+j];
sum += x[i*N+2] * y[2*N+j];
sum += x[i*N+3] * y[3*N+j];
sum += x[i*N+4] * y[4*N+j];
sum += x[i*N+5] * y[5*N+j];
sum += x[i*N+6] * y[6*N+j];
sum += x[i*N+7] * y[7*N+j];
sum += x[i*N+8] * y[8*N+j];
sum += x[i*N+9] * y[9*N+j];
#else
forall (k, N)
sum += x[i*N+k] * y[k*N+j];
#endif
z[i*N+j] = sum;
}
}
当我这样做的时候，我得到了：
VAR=3 USE_STACK=0 N=2000 (unknown) N2=10 (unknown)
Elapsed: 0:00:07.48 (compared with 10.39 / 6.86, not bad, more may be possible).
因此这是分析此类性能问题所需的过程，您需要好的工具。我不知道你的情况发生了什么，因为我无法重现这个问题，但当(小)循环计数未知时，循环展开是MSVC的主要因素(正如预期的那样)。
我使用的测试代码在这里，以防有人想引用它。我想你应该给我一个投票，OP.
编辑：
在Wandbox玩了一段时间，gcc为9.0.0。计时(由于我们在共享盒子上运行，或者更有可能在虚拟机中运行，所以这些计时比较慢，也有点不精确)：
VAR=0 USE_STACK=0 N=2000(已知)N2=10(已知)运行时间=~8秒
VAR=3 USE_STACK=0 N=2000(未知)N2=10(未知)运行时间=约15.5秒
VAR=3 USE_STACK=0 N=2000(未知)N2=10(未知)，循环展开经过时间=~13.5秒
因此，这需要更多的调查——通过探查器和查看生成的代码——而且距离OP的结果还有一百万英里。
现场演示-如果你想尝试不同的东西，你可以自己玩
- 为什么在读取文件大小时文件IO速度会发生变化
- 重载运算符new[]的行为取决于析构函数
- 为什么std::condition_variable notify_all的工作速度比notify_one快(对于随机请
- 文件系统:复制功能的速度秘诀是什么
- 学习多线程C++：添加线程不会使执行速度更快，即使它看起来应该
- 在C++中使用并行化的预期速度是多少(不是 OpenMp，而是 <thread>)
- 两个连续的 OpenMP 并行区域会相互减慢速度
- 新的放置取决于 iostream
- Writefile(） 无法写入数据，具体取决于数据的长度
- 查找标准::hash_map与标准::矢量的速度
- 加快在C++中读取/处理日志文件的速度
- ASIO signal_set多个 IO 线程不可靠，具体取决于代码顺序?
- 为什么这些算法的运行速度比它们应该的要快?
- SFINAE是否取决于类型推断?
- 将强制转换简化为取决于参数的类型
- 修复"-Wunused-parameter"取决于预处理器条件的警告
- 如何提高文件的读取速度?
- 通过libpqxx提高PostgreSQL数据库的更新速度
- 矩阵乘法的速度取决于愚蠢的事情
- 对齐数据:速度取决于恒定数据读取的相对采样大小和采样频率
