Скорость умножения матриц зависит от глупостей

#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

Вступительное примечание: Большая часть этого поста была переписана, поэтому некоторые комментарии ниже не будут иметь особого смысла. Пожалуйста, поищите подробности, если вам это интересно. Так…

ТЛ; др

• профиль, чтобы найти горячие петли
• используйте постоянные значения для них, если это возможно
• попробуйте вручную развернуть их, если нет
• Результаты ОП загадка, никто больше не может воспроизвести что-либо столь экстремальное

Я согласен с @ user4581301 — чем больше компилятор знает заранее, тем больше он может сделать для вас с точки зрения оптимизации вашего кода.

Но вам нужно профилировать этот код — время настенных часов займет у вас столько времени. Я не очень много знаю о профилировщиках для gcc (у ​​меня есть хороший для MSVC), но вы можете попытать счастья Вот.

Это также платит (как @RetiredNinja сказал сразу), чтобы попытаться изучить некоторый ассемблер, используя Godbolt в качестве инструмента, особенно если вы хотите понять столь резкое замедление.

Сказав все это, ваши сроки не имеют смысла для меня, поэтому в вашей системе происходит что-то странное. Поэтому я провел несколько собственных тестов, и мои результаты заметно отличаются от ваших. Я выполнил некоторые из этих тестов на MSVC (потому что там у меня есть такие замечательные инструменты профилирования), а некоторые на gcc на Mac (хотя я думаю, что это на самом деле тайно лязгает под капотом). У меня нет linux box, извините.

Во-первых, давайте разберемся с проблемой размещения таких больших объектов в стеке. Это явно не разумно, и я не могу сделать это вообще на 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;

=> Forall (K, N)
=> сумма + = x [i * N + k] * y [kN + J];
г [яN + 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, когда (небольшое) число циклов неизвестно.

Код теста, который я использовал Вот в случае, если кто-то хочет сослаться на это. Я думаю, ты должен мне проголосовать, ОП.

Редактировать:

Немного поиграл в 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.

Живая демо — вы можете поиграть с этим сами, если хотите попробовать разные вещи, OP.