Question

Я написал программу для быстрого умножения матриц. Чтобы максимально использовать кэш ЦП, он делит матрицу на 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 не играет никакой роли, если мы удалим ее, результаты будут такими же.

Paul Sanders · Answer 1 · 30 июня 2018

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

Т.Л., др

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

Я согласен с @ 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) => sum + = x [i * N + k] * y [k N + j]; z [i 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, когда (небольшое) число циклов неизвестно.

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

Edit:

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

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

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

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

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Похожие темы