Комментарии 44

Интересно было бы сравнить по скорости с библиотекой MKT от того же Intel,
и например с Eigen еще.

На больших матрицах думаю будет небольшой проигрыш в пределах 5-10% (для 1 — потока). На малых размерах проигрыш будет в разы (причины этого указаны в статье).
А померять? Штрассен должен давать ускорение в 10-100 раз на больших матрицах и в 1000+ на ОЧЕНЬ больших и хорошие библиотеки обязаны его поддерживать в теории
UPD: вижу ваш ответ ниже. Странно.
Штрассен на больших матрицах будет накапливать ошибки округления
А почему даже после введения правильного выравнивания используются инструкции для чтения/записи по невыровненному адресу? Разве от них не будет оверхеда? И вообще, мне кажется лучшим решением будет просто ввести требование по выравниванию непосредственно для входных данных.
Я уже достаточно давно разработываю проект по оптимизации различных алгоритмов (в основном по компьютерному зрению) при помощи SIMD инструкций. В частности я всегда старался выравнивать данные и использовать варианты с выровненными инструкциями (например, _mm_load_ps вместо _mm_loadu_ps), где только можно. Однако в последних поколениях процессоров Интел скорость работы этих инструкций одинаковая (при условии обращения к выровненному адресу конечно). Более того компилятор часто вместо _mm_load_ps подставляет ассемблерную инструкцию соответсвующую _mm_loadu_ps.

Быстрое гугление говорит, что разница между _mm_load_p и _mm_loadu_p на выравненных данных исчезла начиная с Nehalem (2008).

Я уже достаточно давно разработываю проект по оптимизации различных алгоритмов (в основном по компьютерному зрению) при помощи SIMD инструкций.


Обёртку на python ещё не завезли? :)
Если подвезете, я возражать не буду. Сам я в этом деле не специалист.
Это была известная отсылка к принципу: «Для всего есть обёртка на python».

Сам я ни разу не писал обёрток, а с графикой работал только в виде лаб в универе, так что ещё больший не специалист, к сожалению.

Очень информативно, а вы пробовали сравнить результаты для реализации "в лоб" у gcc и clang (-O3 -march=haswell -ffast-math)?
Я по правде говоря ожидал что оптимизации из первого шага и так должны быть сделаны компилятором (loop unswitch/rotate, const propogation?)

Вы не пробовали провести что-то подобное для матриц в столбцовом формате?
Основной принцип останется прежним: микроядра, максимальная локализация в кэше. Изменится только функции переупорядочивания.
GPU просто созданы для быстрого перемножения матриц.
Зачем пытаться их обогнать на CPU? Ведь все равно не выйдет.
Цель статьи — не написать саму быструю реализацию, а показать внутренее устройство алгоритма. Понятно, что GPU будет быстрее. С этим никто не спорит.
А где можно почитать (или может, у вас есть в планах написать статью) об особенностях реализации умножения на картах?

Интересно, какие там аппаратные особенности и как они используются.

Спасибо.
Я бы упомянул симпатичную математическую абстракцию часто используемую в подобных в подобных оптимизациях. Для матричного умножения справедливо что если мы поделим матрицу на блоки меньшего размера, то матрицы можно перемножать перемножая блоки

Это можно использовать либо явно (tiled версия aka блочное перемножение матриц) либо в рекурсивных схемах (где возникают всякие вкусные субкубические варианты типа алгоритма Штрассена). У Вас как я понимаю по сути tiled версия сильно заточенная под AVX — ширина блоков выбирается под ширину регистров, длина — под размер L1.
На сколько я видел в разных реализациях, нигде блочное перемножение не используется. Видимо потому, что это приводит к дополнительным копированиям памяти. А как следует из статьи выигрыш от локализации данных в кэше — до 10 раз. Выигрыш от подобных блочных алгоритмов будет проявляться только на матрицах циклопических размеров. А такие уже будет целесообразно считать на GPU.
На платформе РС может быть и не используется, а мы в обработке сигналов используем. Хотя ваше микроядро по-моему туда же относится.
Наверное не так выразился. Мысль которую я хотел донести: Не используются субкубические варианты типа алгоритма Штрассена.
Причём важна форма блоков! Типичная задача — перемножить большие матрицы в медленной памяти, цикл такой: разбиваем на блоки, пересылаем блоки аргументов во внутреннюю память, перемножаем, выгружаем блок результата. Как минимизировать пересылки? Оказывается, чем крупнее и квадратнее блок, тем лучше, поскольку объём работы по перемножению матриц IxK и KxJ пропорционален I*J*K, а объём пересылок — I*K + K*J + I*J.
Вспомнилось, как я году эдак в 2013 матрицы перемножал. Правда, там в ячейках были не float, а int64 (точнее, в A*B=C в A и B — инты, а в C должно было инт64 влезать). Размеры матричек 5000х5000. Подробнее о задаче на fastcomputing.org (вебархив, ибо проект ныне покойный). А вот табличка рекордов. В 40 секунд тогда упихал. Подробности реализации сейчас не вспомню, но вроде там Виноград-Штрассен, который где-то то ли на 200х200, то ли на 1000х1000 переключается уже на перемножение в лоб. Без Штрассена было тормознее, насколько помню. В перемножении в лоб тоже всякие SSE, помню разве что только один хак: мы умножаем не строчку A на строчку B (B транспонирована, конечно), а пару строчек A на пару строчек B, тем самым сокращая число чтений, потом все идет чисто на регистрах процессора, и на выходе получаем 4 значения для этих самых пар строчек.
MKL конечно платный, не всем подходит, но есть же MKL-DNN открытый.
Для мелких лучше использовать libxsmm. Он делает JIT компиляцию микрокернела под вашу платформу и бесплатный.

Тем более. Берите МКЛ и забудьте про перемножение матриц. По крайней мере кроме краевых примеров 2х2.
Наверное, стоит добавить, что когда матриц много, то очень часто имеет смысл сначала определить порядок умножения матриц.
Очень полезная и хорошая статья. Насколько я знаю, MKL всегда работает на числе потоков, равном числу физических ядер. Насколько результаты статьи зависят от числа использованных виртуальных/физических потоков?
Цитата из введения: С целью ограничить объем изложения, я ограничился описанием однопоточного алгоритма для обычных процессоров. Тема многопоточности и алгоритмов для графических ускорителей явно заслуживает отдельной статьи.
Уже увидел. Значит это вопрос к следующей статье. С MKL, как с универсальным средством сравнения для высокопроизводительных вычислений, полезно приводить сравнение всегда.
MKL работает:
а) Если выставлен MKL_NUM_THREADS — то количество потоков ограничено этой переменной;
б) Если выставлен OMP_NUM_THREADS и не выставлена а) то ограничивается этой переменной;
По умолчанию считается что количество потоков = количеству ядер. Т.е если ОБЕ переменные не выставлены — будет работать на всем до чего дотянется.
Т.е. для использования всех виртуальных ядер с учетом гипертрейдинга достаточно установить MKL_NUM_THREADS или OMP_NUM_THREADS в удвоенное количество физических ядер? Заранее спасибо.
С точки зрения ОС — гипертрединг это и так двойное кол-во ядер, так что вообще не ставьте переменную, оно само определит.
Может быть и определяет, но использует только физические ядра.
Попробовал на dgemm с включенным гипертрейдингом. Уменьшить число потоков через set MKL_NUM_THREADS получается, а увеличить больше, чем число физических ядер — нет. Быстродействие упрямо показывает мах 50% загруженности cpu. Аналогично mkl_set_num_threads не дает загрузить число ядер больше физического. Функция mkl_get_max_threads возвращает значение равное числу физических ядер, а не виртуальных. Помнится, разработчики из Нижнего Новгорода подтверждали этот факт.
Тут решение уже за MKL. То, что вы видите загрузку в 50% — к реальности может(и скорее всего так и есть) не иметь никакого отношения, так как каждое из виртуальных ядер может полностью загрузить реальное. Т.е. с точки зрения операционки, одно ядро загружено полностью, на втором вообще ничего не запущено — загрузка 50%, а на деле единственное физическое ядро имеет полностью загруженные конвееры FPU.
MKL библиотека от производителя процессора, я сильно сомневаюсь чтобы они специально занижали производительность своих продуктов.
А в случае когда два потока работают на одном физическом ядре они разделяют как блоки FPU, так и кеш, что скажется на производительности скорее негативно, чем позитивно.
Использование гипертридинга в нагруженных кернелах не приведет ни к чему хорошему.
50% загрузка видимо видна в диспетчере процессов Windows, там да. Логическое ядро отображается как нормальное, в связи с чем получается такое недопонимание.
Запустив одновременно две задачи, убедился, что время расчета выросло в два раза.
Статья читается как какая-то магия. Я, как джавист из энтерпрайза, смог бы самостоятельно придти только к первому варианту решения, и если бы после этого встал вопрос повышения производительности, то только и смог бы придумать что «попробовать в несколько потоков выполнять» или «раскидать вычисления по облаку».
Делал подобное лет 10 назад, использовал практически все те же оптимизации, правда оптимизацией под L2/L3 уже не заморачивался, прирост не особо большой от этого получался. Результатами до сих пор пользуюсь — библиотека получилась небольшой, не надо что-то стороннее тянуть большое.
Но все-таки умножение матриц это достаточно простая вещь. Вот сделать то же самое для обращения матрицы (LUP разложения) и распараллелить обращение было посложнее.
А что если перед умножением транспонировать матрицу B? Тогда по идее можно получить ускорение за счет того, что элементы матрицы B, которые участвуют в умножении на одну строку матрицы A станут соседями в памяти и будут лучше размещаться в кеше.
К сожалению нет. На первый взгляд действительно все проще: данные A и B лежат одинаково — только считай взаимное скалярное произведение их строчек.

Однако: максимальный размер микроядра получится 3x4, что дает нам (3 + 4)/(3*4) = ~0.58 загрузок на одну fma. Напомню, что при классической схеме с окном 6x16 получается (6 + 16)/(6*16) = ~0.23 загрузок на одну fma. Т.е. предложенная вами схема почти в 2.5 раза более требовательна к пропускной способности памяти. В принципе мои внутренние тесты это подтверждают.
Спасибо, классная статья. Не ожидал, что можно так сильно оптимизировать.

В начале статьи Вы критикуете малопонятность опенсурсных реализаций, ссылаясь в большей степени на ассемблер. Могу покритиковать Ваш код, так как он использует захардкоженные значения :) Думаю, именованные константы мне бы были более понятны, чем всякие 6, 24 и т. д…
Как раз имею Core i5 4-го поколения (должно иметь приличные FMA), правда использовал для Matlab и для FFT (вроде как там по такой же логике польза от FMA) на комплексных числах. Вроде как прирост 4670 против Cel. 1005M достаточно большой был.
Ещё раньше использовал Athlon II X3 445, но для немного других задач. Там можно было явно записать решение задачи
dE/dt = f(E)*E + a*(laplas)E
в форме
dE/E = [f(E) — a*k2] dt.
То есть это обычный метод с функцией fft. При этом 3-ядерник грузился процентов на 60, но брат думает, что это тратилось второе ядро только на вывод графики.
Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.