Почему наивное умножение матриц C++ в 100 раз медленнее, чем BLAS?

Я изучаю умножение больших матриц и провел следующий эксперимент, чтобы сформировать базовый тест:

1. Произвольно сгенерируйте две матрицы 4096x4096 X, Y из стандартного нормального (0 среднее значение, 1 стандартное отклонение).
2. Z = X*Y
3. Суммируйте элементы Z (чтобы убедиться, что они доступны) и выведите.

Вот наивная реализация 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

Octave под капотом использует BLAS (я предполагаю функцию sgemm)

Аппаратное обеспечение — i7 3930X на Linux x86-64 с 24 ГБ оперативной памяти. BLAS использует два ядра. Возможно, гиперпоточная пара?

Я обнаружил, что версия C++, скомпилированная с помощью GCC 4.7 в -O3, выполнялась 9 минут:

real    9m2.126s
user    9m0.302s
sys         0m0.052s

Октавная версия заняла 6 секунд:

real    0m5.985s
user    0m10.881s
sys         0m0.144s

Я понимаю, что BLAS оптимизирован до чертиков, а наивный алгоритм полностью игнорирует кеши и прочее, но серьезно -- 90 раз?

Кто-нибудь может объяснить эту разницу? Какова именно архитектура реализации BLAS? Я вижу, что он использует Фортран, но что происходит на уровне процессора? Какой алгоритм он использует? Как он использует кеши процессора? Какие машинные инструкции x86-64 он вызывает? (Использует ли он расширенные функции ЦП, такие как AVX?) Откуда он получает эту дополнительную скорость?

Какие ключевые оптимизации алгоритма 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% наблюдаемое использование ЦП или нет.

Какой библиотекой является 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

Это АТЛАС 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)