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

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

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

Вот наивная реализация С++:

#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))

Run:

$ time octave < test.octave

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

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

Я обнаружил, что версия С++, скомпилированная с 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? Я вижу, что он использует Fortran, но что происходит на уровне CPU? Какой алгоритм он использует? Как это работает с кэшами процессора? Что это за машинные инструкции x86-64? (Использует ли он расширенные функции процессора, такие как AVX?) Откуда у вас эта дополнительная скорость?

Какие ключевые оптимизации для алгоритма С++ могут получить его на одном уровне с версией 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

Это ATLAS 3.8.4

Вот оптимизаторы, которые я позже реализовал:

Использование элемента с черепицей, где я предварительно загружаю 64x64 блоков X, Y и Z в отдельные массивы.

Изменение вычисления каждого блока так, чтобы внутренний цикл выглядел так:

for (size_t tcol = 0; tcol < block_width; tcol++)
    bufz[trow][tcol] += B * bufy[tk][tcol];

Это позволяет GCC автоматически деактивировать инструкции SIMD, а также позволяет уровень команд parallelism (я думаю).

Включение 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)

Ответ 1

Вот три фактора, ответственные за разницу в производительности между вашим кодом и BLAS (плюс примечание к алгоритму Strassens).

В вашем внутреннем цикле, на k, у вас есть y[k*dim + col]. Из-за того, как устроен кеш памяти, последовательные значения k с теми же картами dim и col соответствуют одному и тому же набору кешей. Способ кэширования структурирован, каждый адрес памяти имеет один кеш-набор, где его содержимое должно храниться, пока оно находится в кеше. Каждый набор кешей имеет несколько строк (четыре - это типичное число), и каждая из этих строк может содержать любой из адресов памяти, которые соответствуют этому конкретному набору кешей.

Так как ваш внутренний цикл проходит через y таким образом, каждый раз, когда он использует элемент из y, он должен загружать память для этого элемента в тот же набор, что и предыдущая итерация. Это заставляет одну из предыдущих строк кеша в наборе выходить. Затем в следующей итерации цикла col все элементы y были выселены из кеша, поэтому их необходимо перезагрузить снова.

Таким образом, каждый, когда ваш цикл загружает элемент y, он должен быть загружен из памяти, которая занимает много циклов процессора.

Высокопроизводительный код избегает этого двумя способами. Во-первых, он делит работу на более мелкие блоки. Строки и столбцы разбиты на меньшие размеры и обрабатываются короткими циклами, которые могут использовать все элементы в строке кэша и использовать каждый элемент несколько раз, прежде чем они перейдут к следующему блоку. Таким образом, большинство ссылок на элементы x и элементы y исходят из кеша, часто в одном процессорном цикле. Во-вторых, в некоторых ситуациях код будет копировать данные из столбца матрицы (которая перетаскивает кеш из-за геометрии) в строку временного буфера (что позволяет избежать перекоса). Это снова позволяет использовать большинство ссылок на память из кеша вместо памяти.

Другим фактором является использование функций Single Instruction Multiple Data (SIMD). Многие современные процессоры имеют инструкции, которые загружают несколько элементов (четыре элемента float являются типичными, но некоторые из них теперь выполняют восемь) в одной инструкции, хранят несколько элементов, добавляют несколько элементов (например, для каждого из этих четырех, добавляют к соответствующему из этих четырех), умножить несколько элементов и т.д. Простое использование таких инструкций сразу же делает ваш код в четыре раза быстрее, если вы можете организовать свою работу для использования этих инструкций.

Эти инструкции напрямую не доступны в стандартном C. Некоторые оптимизаторы теперь пытаются использовать такие инструкции, когда могут, но эта оптимизация затруднена, и нередко получается, что она очень выигрывает от нее. Многие компиляторы предоставляют расширения для языка, предоставляющего доступ к этим инструкциям. Лично я обычно предпочитаю писать на ассемблере, чтобы использовать SIMD.

Другим фактором является использование функций параллельного выполнения на уровне инструкций на процессоре. Обратите внимание, что в вашем внутреннем цикле обновляется acc. Следующая итерация не может добавить к acc до тех пор, пока предыдущая итерация не завершит обновление acc. Высокопроизводительный код вместо этого будет содержать несколько сумм, работающих параллельно (даже несколько SIMD-сумм). Результатом этого будет то, что, пока выполняется добавление для одной суммы, будет добавлено добавление для другой суммы. На современных процессорах принято поддерживать четыре или более операций с плавающей запятой одновременно. Как написано, ваш код не может этого сделать вообще. Некоторые компиляторы будут пытаться оптимизировать код, переставляя петли, но это требует, чтобы компилятор мог видеть, что итерации определенного цикла независимы друг от друга или могут быть заменены другим циклом и т.д.

Вполне возможно, что использование кеша эффективно обеспечивает десятикратное повышение производительности, SIMD предоставляет еще четыре, а уровень parallelism на уровне инструкций предоставляет еще четыре, давая 160 в целом.

Вот очень грубая оценка эффекта алгоритма Страссенса на основе этой страницы Википедии. На странице Википедии Штрассен немного лучше прямого умножения вокруг n = 100. Это говорит о том, что отношение постоянных факторов времени выполнения составляет 100 ³/100 ^2.807 ≈ 2.4, Очевидно, что это будет сильно варьироваться в зависимости от модели процессора, размеров матриц, взаимодействующих с эффектами кеша, и так далее. Однако простая экстраполяция показывает, что Штрассен примерно в два раза эффективнее прямого умножения при 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 linkage (по умолчанию), тогда компилятор должен разрешить которая вызывается из внешнего модуля, поэтому она должна генерировать код, который работает правильно, если bufz и bufy перекрываются. (Один из способов, с которым компилятор может справиться с этим, состоит в том, чтобы сгенерировать две версии подпрограммы, одну из которых вызывать из внешних модулей, и одну из них, которая будет вызываться из модуля, который в настоящее время компилируется. Независимо от того, зависит ли это от вашего компилятора, и т.д.). Если вы объявляете подпрограмму как static, то компилятор знает, что ее нельзя вызывать из внешнего модуля (если вы не берете его адрес, и есть возможность, что адрес передается вне текущего модуля).

Еще одна потенциальная проблема заключается в том, что даже если компилятор векторизует этот код, он не обязательно генерирует лучший код для процессора, который вы выполняете. Рассматривая сгенерированный код сборки, кажется, что компилятор использует только %ymm1 несколько раз. Снова и снова она умножает значение из памяти на %ymm1, добавляет значение из памяти в %ymm1 и сохраняет значение из %ymm1 в память. Есть две проблемы с этим.

Во-первых, вы не хотите, чтобы эти частичные суммы сохранялись в памяти часто. Вы хотите, чтобы все дополнения были записаны в регистр, и регистр будет записан в память нечасто. Убеждение компилятора в этом, вероятно, требует переписывания кода для явного указания частичных сумм во временных объектах и записи их в память после завершения цикла.

Два, эти инструкции номинально зависят последовательно. Добавление не может начинаться до тех пор, пока умножение не завершится, и хранилище не сможет записать в память до завершения добавления. Core i7 обладает большими возможностями для исполнения вне порядка. Таким образом, хотя у этого есть, что добавляет ожидание, чтобы начать выполнение, он смотрит на умножение позже в потоке команд и запускает его. (Даже при том, что это умножение также использует %ymm1, процессор перенаправляет регистры "на лету", так что для этого используется другой внутренний регистр.) Несмотря на то, что ваш код заполнен последовательными зависимостями, процессор пытается выполнить несколько инструкции сразу. Однако некоторые вещи могут помешать этому. Вы можете исчерпать внутренние регистры, которые процессор использует для переименования. Используемые вами адреса памяти могут приводить к ложным конфликтам. (Процессор просматривает около десятка бит младших бит адресов памяти, чтобы узнать, может ли адрес быть таким же, как другой, который он пытается загрузить или сохранить из более ранней инструкции. Если бит равен, процессор имеет чтобы задержать текущую нагрузку или сохранить до тех пор, пока она не сможет проверить, что весь адрес отличается. Эта задержка может увеличиваться больше, чем только текущая загрузка или сохранение.) Таким образом, лучше иметь инструкции, которые явно независимы.

Это еще одна причина, по которой я предпочитаю писать высокопроизводительный код в сборке. Чтобы сделать это на C, вам нужно убедить компилятор дать вам такие инструкции, например, написать некоторые из ваших собственных SIMD-кодов (с использованием расширений языка для них) и вручную развернуть петли (выписывая несколько итераций).

При копировании в буфер и из него могут возникать аналогичные проблемы. Тем не менее, вы сообщаете, что 90% времени потрачено на calc_block, поэтому я не смотрел на это внимательно.

Кроме того, алгоритм Strassens учитывает любую оставшуюся разницу?

Ответ 2

Алгоритм Strassen имеет два преимущества по сравнению с наивным алгоритмом:

Лучшая временная сложность с точки зрения количества операций, как правильно указывают другие ответы;
Это кэш-забывающий алгоритм. Разница в количестве промахов кэша находится в порядке B*M½, где B - размер строки кэша, а M - размер кэша.

Я думаю, что второй пункт объясняет многое для замедления, которое вы испытываете. Если вы запускаете свое приложение под Linux, я предлагаю вам запустить их с помощью инструмента perf, в котором указывается, сколько промахов в кеше выполняется программой.

Ответ 3

Я не знаю, насколько достоверна информация, но Wikipedia говорит, что BLAS использует алгоритм Strassen для больших матриц. И твой большой. Это около O (n ^ 2,807), что лучше, чем ваш наивный алорифм O (n ^ 3).

Ответ 4

Это довольно сложная тема, и на этот вопрос ответил Эрик. Я просто хочу указать на полезную ссылку в этом направлении, стр. 84:

http://www.rrze.fau.de/dienste/arbeiten-rechnen/hpc/HPC4SE/

который предлагает сделать "цикл разворачивания и застревания" поверх блокировки.

Может ли кто-нибудь объяснить эту разницу?

Общее объяснение состоит в том, что отношение числа операций/числа данных равно O (N ^ 3)/O (N ^ 2). Таким образом, матрично-матричное умножение является алгоритмом, связанным с кешем, что означает, что вы не страдаете общим узким местом пропускной способности для больших размеров матрицы. Если код хорошо оптимизирован, вы можете получить до 90% максимальной производительности вашего ЦП. Таким образом, потенциал оптимизации, разработанный Эриком, огромен, как вы заметили. На самом деле было бы очень интересно увидеть лучший исполняемый код и скомпилировать вашу финальную программу с другим компилятором (обычно похвастаться интеллектом).

Ответ 5

Около половины разницы учитывается при алгоритмическом улучшении. (4096 * 4096) ^ 3 - сложность вашего алгоритма, или 4.7x10 ^ 21, и (4096 * 4096) ^ 2.807 - 1x10 ^ 20. Это разница около 47 раз.

Другие 2x будут учитываться более интеллектуальным использованием кеша, инструкций SSE и других подобных низкоуровневых оптимизаций.

Изменить: я ложь, n - ширина, а не ширина ^ 2. Алгоритм на самом деле будет составлять около 4x, так что еще около 22x идти. Инструкции потоков, кеша и SSE могут хорошо учитывать такие вещи.