Почему умножение наивной матрицы C ++ в 100 раз медленнее, чем BLAS?
Я взглянул на умножение больших матриц и провел следующий эксперимент, чтобы сформировать базовый тест:
- Произвольно сгенерируйте две матрицы 40, x 4096 X, Y из стандартного стандартного значения (среднее значение 0, стандартное значение 1).
- Z = X * Y
- Суммируйте элементы 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
Октав под капотом использует 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? Я вижу, что он использует Fortran, но что происходит на уровне процессора? Какой алгоритм он использует? Как это с использованием кешей процессора? Какие машинные инструкции для 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
Это ATLAS 3.8.4
Вот оптимизации, которые я позже реализовал:
Используя плиточный подход, где я предварительно загружаю 64×64 блоки X, Y и Z в отдельные массивы.
Изменение расчета каждого блока так, чтобы внутренний цикл выглядел так:
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)
Решение
Вот три фактора, ответственных за разницу в производительности между вашим кодом и BLAS (плюс примечание об алгоритме Штрассена).
В вашем внутреннем цикле, на k, у тебя есть y[k*dim + col], Из-за того, как устроен кэш памяти, последовательные значения k с тем же dim а также col сопоставить с тем же набором кэша. Как устроен кэш, каждый адрес памяти имеет один набор кеша, в котором его содержимое должно храниться, пока он находится в кеше. Каждый набор кеша имеет несколько строк (четыре — типичное число), и каждая из этих строк может содержать любой из адресов памяти, которые сопоставляются с этим конкретным набором кеша.
Потому что ваш внутренний цикл проходит через y таким образом, каждый раз, когда он использует элемент из y, он должен загрузить память для этого элемента в тот же набор, что и предыдущая итерация. Это заставляет удалить одну из предыдущих строк кэша в наборе. Затем в следующей итерации col цикл, все элементы y были удалены из кэша, поэтому они должны быть перезагружены снова.
Таким образом, каждый время ваш цикл загружает элемент y, он должен быть загружен из памяти, что занимает много тактов процессора.
Высокопроизводительный код позволяет избежать этого двумя способами. Во-первых, он делит работу на более мелкие блоки. Строки и столбцы разделены на меньшие размеры и обработаны с помощью более коротких циклов, которые могут использовать все элементы в строке кэша и использовать каждый элемент несколько раз, прежде чем перейти к следующему блоку. Таким образом, большинство ссылок на элементы x и элементы y приходят из кеша, часто в одном цикле процессора. Во-вторых, в некоторых ситуациях код будет копировать данные из столбца матрицы (которая разбивает кэш из-за геометрии) в строку временного буфера (который избегает перебора). Это опять-таки позволяет использовать большинство ссылок на память из кеша, а не из памяти.
Другим фактором является использование функций единой инструкции нескольких данных (SIMD). Многие современные процессоры имеют инструкции, которые загружают несколько элементов (четыре float элементы типичны, но некоторые теперь выполняют восемь) в одной инструкции, сохраняют несколько элементов, добавляют несколько элементов (например, для каждого из этих четырех, добавляют его к соответствующему одному из этих четырех), умножают несколько элементов и так далее. Простое использование таких инструкций сразу делает ваш код в четыре раза быстрее, если вы можете организовать свою работу для использования этих инструкций.
Эти инструкции не доступны напрямую в стандарте C. Некоторые оптимизаторы сейчас пытаются использовать такие инструкции, когда могут, но эта оптимизация трудна, и нечасто получить от нее большую выгоду. Многие компиляторы предоставляют расширения для языка, которые дают доступ к этим инструкциям. Лично я обычно предпочитаю писать на ассемблере, чтобы использовать SIMD.
Другим фактором является использование функций параллельного выполнения на уровне команд на процессоре. Обратите внимание, что в вашем внутреннем цикле, acc обновляется. Следующая итерация не может добавить к acc пока предыдущая итерация не закончила обновление acc, Высокопроизводительный код вместо этого будет поддерживать одновременную работу нескольких сумм (даже нескольких сумм SIMD). Результатом этого будет то, что во время выполнения сложения для одной суммы будет начато сложение для другой суммы. На современных процессорах принято поддерживать одновременно четыре или более операций с плавающей запятой. Как написано, ваш код не может сделать это вообще. Некоторые компиляторы пытаются оптимизировать код, переставляя циклы, но для этого необходимо, чтобы компилятор мог видеть, что итерации определенного цикла независимы друг от друга или могут быть заменены другим циклом и так далее.
Вполне возможно, что эффективное использование кэша обеспечивает повышение производительности в 10 раз, SIMD — еще четыре, а параллелизм на уровне команд — еще четыре, что в сумме дает 160.
Вот очень грубая оценка эффекта алгоритма Штрассена, основанная на эта страница Википедии. На странице Википедии сказано, что Штрассен немного лучше, чем прямое умножение вокруг n = 100. Это говорит о том, что отношение постоянных факторов времени выполнения равно 1003 / 1002,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 связи (по умолчанию), то компилятор должен разрешить вызов подпрограммы из внешнего модуля, поэтому он должен генерировать код, который работает правильно, если bufz а также bufy перекрытия. (Один из способов, которым компилятор может справиться с этим, — это сгенерировать две версии подпрограммы, одну для вызова из внешних модулей и одну для вызова из модуля, который компилируется в настоящее время. Зависит ли это от вашего компилятора, переключателей оптимизации, и так далее.) Если вы объявите процедуру как staticтогда компилятор знает, что он не может быть вызван из внешнего модуля (если вы не берете его адрес и есть вероятность, что адрес будет передан за пределы текущего модуля).
Другая потенциальная проблема заключается в том, что даже если компилятор векторизует этот код, он не обязательно генерирует лучший код для процессора, на котором вы выполняете. Глядя на сгенерированный ассемблерный код, кажется, что компилятор использует только %ymm1 несколько раз. Снова и снова, он умножает значение из памяти в %ymm1, добавляет значение из памяти в %ymm1и сохраняет значение из %ymm1 в память. Есть две проблемы с этим.
Во-первых, вы не хотите, чтобы эти частичные суммы часто сохранялись в памяти. Вы хотите, чтобы в регистр накапливалось много дополнений, и регистр будет записываться в память нечасто. Чтобы убедить компилятор сделать это, вероятно, потребуется переписать код, чтобы он явно содержал частичные суммы во временных объектах и записывал их в память после завершения цикла.
Во-вторых, эти инструкции являются номинально серийно зависимыми. Добавление не может начаться до тех пор, пока умножение не завершится, и хранилище не сможет записать в память, пока не завершится добавление. Core i7 обладает отличными возможностями для выполнения заказов вне очереди. Таким образом, пока у него есть это добавление, ожидающее начала выполнения, оно просматривает умножение позже в потоке инструкций и запускает его. (Хотя это умножение также использует %ymm1процессор перераспределяет регистры на лету, так что он использует другой внутренний регистр для этого умножения.) Даже если ваш код заполнен последовательными зависимостями, процессор пытается выполнить несколько инструкций одновременно. Тем не менее, некоторые вещи могут помешать этому. Вы можете исчерпать внутренние регистры, которые процессор использует для переименования. Адреса памяти, которые вы используете, могут столкнуться с ложными конфликтами. (Процессор просматривает дюжину или около того младших битов адресов памяти, чтобы увидеть, может ли адрес быть таким же, как другой, который он пытается загрузить или сохранить из более ранней инструкции. Если биты равны, процессор имеет задержать текущую загрузку или сохранение до тех пор, пока не сможет проверить, что весь адрес отличается. Эта задержка может запутать больше, чем просто текущая загрузка или сохранение.) Поэтому лучше иметь инструкции, которые явно независимы.
Это еще одна причина, по которой я предпочитаю писать высокопроизводительный код на ассемблере. Чтобы сделать это в C, вы должны убедить компилятор дать вам такие инструкции, выполнив такие вещи, как написание некоторого собственного кода SIMD (используя для них языковые расширения) и ручное развертывание циклов (запись нескольких итераций).
При копировании в буферы и из них могут возникнуть похожие проблемы. Тем не менее, вы сообщаете, что 90% времени тратится на calc_block, поэтому я не смотрел на это внимательно.
Кроме того, учитывает ли алгоритм Штрассена какое-либо из оставшихся отличий?
Другие решения
Алгоритм Штрассена имеет два преимущества по сравнению с наивным алгоритмом:
- Лучшая временная сложность с точки зрения количества операций, как правильно указывают другие ответы;
- Это алгоритм кеширования. Разница в количестве промахов кеша в порядке
B*M½где B — размер строки кэша, а M — размер кэша.
Я думаю, что второй момент объясняет замедление, которое вы испытываете. Если вы запускаете свое приложение под Linux, я предлагаю вам запустить их с perf инструмент, который говорит вам, сколько пропусков кэша программа испытывает.
Я не знаю, насколько достоверна информация, но Википедия говорит, что BLAS использует алгоритм Штрассена для больших матриц. И ваши действительно большие. Это примерно O (n ^ 2.807), что лучше, чем ваш O (n ^ 3) наивный алогритм.
Это довольно сложная тема, и Эрик хорошо ответил в посте выше. Я просто хочу указать на полезную ссылку в этом направлении, стр. 84:
http://www.rrze.fau.de/dienste/arbeiten-rechnen/hpc/HPC4SE/
который предлагает сделать «развернуть петлю и замять» поверх блокировки.
Кто-нибудь может объяснить эту разницу?
Общее объяснение состоит в том, что отношение числа операций / количества данных составляет O (N ^ 3) / O (N ^ 2). Таким образом, умножение матрицы на матрицу является алгоритмом с привязкой к кешу, что означает, что вы не страдаете от общего узкого места в пропускной способности памяти для больших размеров матриц.
Вы можете получить до 90% максимальной производительности вашего процессора, если код хорошо оптимизирован. Так что потенциал оптимизации, разработанный Эриком, огромен, как вы заметили. На самом деле, было бы очень интересно увидеть наиболее эффективный код и скомпилировать финальную программу с другим компилятором (Intel обычно хвастается, что он лучший).
Около половины различий приходится на улучшение алгоритмов. (4096 * 4096) ^ 3 — сложность вашего алгоритма, или 4,7×10 ^ 21, а (4096 * 4096) ^ 2.807 — 1×10 ^ 20. Это разница примерно в 47 раз.
Другие 2x будут объяснены более разумным использованием кэша, инструкций SSE и других подобных низкоуровневых оптимизаций.
Редактировать: я лгу, n это ширина, а не ширина ^ 2. Алгоритм на самом деле будет составлять только около 4x, так что еще предстоит 22x. Потоки, кэш и инструкции SSE вполне могут объяснить такие вещи.