Быстрое перемножение матриц

Автор: Zealint, alexBlack - 22.02.2010
Редактор: Максим Лович
Оригинал Zealint. Оригинал alexBlack.

Товарищ Zealint провел конкурс по перемножению матриц, в котором нужно было перемножить две матрицы размером 5000×5000.

В конкурсе приняли участие 8 человек. С большим отрывом от остальных победил тов. venco.

1. venco :: C++ GCC 4.0.2 :: 7,52 с.
2. blackmirror :: fasm :: 11,13 с.
3. Zealint :: Intel C++ 11.1 + ASM (Виноград — Штрассен) :: 11,33 с.
4. alexBlack :: Delphi (ручн. реал.) + MMX :: 19,2 с.
5. ghoul :: GCC 3.4.5 + gotoBLAS2 1.13 :: 35,84 с.
6. Ulex :: Masm 6.14.8444 :: 57,74 с.
7. covax :: MS VC 2010 :: 62,89 с.
8. yakasuka :: VS6 (ручн. реал.) :: 144,44 с.

Решение тов. Zealint (исходники в архиве)
(цитата)

Во-первых, сначала я написал самую тупую версию, которая умножает строку на столбец в тройном цикле (C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ], циклы идут в порядке i-j-k). Это и были те самые 800 с. (надо же с чего-то начинать). Здесь же можно поиграть перестановкой циклов местами: самое быстрое получается, если выбрать i-k-j, в этом случае все три матрицы сканируются по строкам. Такая версия уже укладывалась в 3 минуты (что почти в 5 раз быстрее). Можно просто транспонировать вторую матрицу прямо на вводе и поменять местами индексы k и j при обращении к матрице B. Это ещё быстрее (150 с.). Что ещё можно сделать без вмешательства ассемблера? — поменять int (32 бита) на short (16 бит) для входящих матриц и скомпилировать родным для процессора компилятором с опциями максимального ускорения. Это дает 94 секунды. До 55 секунд можно дойти, если развернуть матрицы в массив. Тогда будут получаться записи вида A [ i * n + k ], но все такие умножения нужно либо вынести за внутренний цикл, либо вообще заменить сложением указателей с числом n в нужных местах. В этом случае компилятор оптимизирует гораздо лучше. Все, без ассемблера можно сделать только ещё одну вещь: разбить матрицу на блоки, которые помещаются в кэш, но я решил, что сделаю это сразу на ассемблере.

Во-вторых, вспоминаем про MMX, там есть инструкция pwaddwd, которая умеет 4 слова умножать на 4 слова и сохранять суммы первых двух умножений и следующих двух. Все это работает на 64-битных mm регистрах и довольно быстро. Не забываем, что размер конвейера — 64 байта, и разворачиваем циклы нужное количество раз. Кстати, я разворачивал циклы так, что получалось больше 64 байт, это совершенно не влияло на время работы, перемешивание инструкций в правильном порядке тоже ничего не даёт, так как процессор сам решает, в каком порядке их исполнять. Все это (вместе с разбиением матрицы на блоки) позволяет получить 33 секунды.

В-третьих, я узнал, что на xmm регистрах операция pwaddwd умеет умножать 8 слов на 8 слов и сохранять 4 двойных слова. Это означает, что за один такт можно перемножить 2×2 матрицу на 2×2 матрицу и сохранить тут же 2×2 матрицу из двойных слов. Итак, блоки, которые помещаются в кэш, нужно тоже перемножать не поэлементно, а блоками 2×2, а для этого нужно правильно считать матрицу. Например, на входе была матрица
1 2 3 4
5 6 7 8
9 10 11 12
13 14 15 16
Нужно считать её так:
1 2 5 6
3 4 7 8
9 10 13 14
11 12 15 16
В один xmm регистр считываем сразу 1, 2, 5, 6, 3, 4, 7, 8, а во второй – то же самое из второй матрицы. То есть мы считали сразу две 2×2 матрицы, а нужно было одну. Перемешиваем элементы так, чтобы было xmm0 = 1 2 1 2 3 4 3 4 и xmm1 = 5 6 5 6 7 8 7 8, аналогично, правильный порядок элементов нужно сделать для xmm2 и xmm3, где хранятся элементы второй матрицы. Теперь за 2 операции pwaddwd мы получаем две матрицы 2×2. Далее, исходные две матрицы 2×2 можно, после умножения на первую строку матрицы B, сразу же умножить на вторую строку матрицы B, таким образом (считайте сами) у меня получилось 0.75 обращений к памяти на каждые 8 умножений (а было — 2 и сначала даже 4). Все это дало 18 секунд. Да, чуть не забыл, массивы нужно выровнять по 16 байт, чтобы использовать инструкцию movdqa. У меня компилятор выравнивал их сам.

В-четвёртых, я изменил скорость ввода и вывода данных. Раньше пользовался gets, теперь решил использовать API функцию ReadFile, это сэкономило еще 4 секунды и получилось 14,5. Это самое лучшее, что удалось выжать из кубической версии.

В-пятых, пишем алгоритм Винограда — Штрассена (7 умножений и 15 сложений) и… у меня не получилось написать его достаточно аккуратно, поэтому только 11,3 секунды. После этого я убрал половину всех вспомогательных матриц, переписал все сложения через SSE, и это совершенно ничего не изменило. Даже чуть медленнее стало. Но до конца конкурса оставалось около часа, и я так и не успел разобраться до конца.

В любом случае меня сильно удивило, что удалось так далеко продвинуться, я не верил, что можно так радикально ускорить программу, полагаясь до этого лишь на компилятор.

… Дальше идут мои собственные соображения по поводу конкурса, ограничений и самой задачи. Все эти размышления основаны на моём личном опыте, поэтому читайте осторожно (не примеряйте на себя и не обижайтесь, если кого задел).

Решение тов. venco
(цитата)

1. Ручной ввод/вывод с использованием отображения файлов в память.
2. Транспорование матрицы B и блочное умножение по 32 строки матрицы B за проход. Это улучшает эффективность кеша. Транспонирование тоже блочное по 128х128 элементов опять же для кеша.
3. Умножение в xmm регистрах инструкцией pmaddwd. Причём я умножаю 4 строки за раз, и результат (4 значения) получаю в одном xmm регистре и сохраняю на выход за раз. Всё выше дало 13 секунд.
4. Штрассен. Сначала на один уровень, а потом до 5-ти. При этом обнаружилось, что сложный вариант с минимумом сложений увеличивает модуль операндов до 4-х раз, что ограничивает глубину 2-мя уровнями, иначе возникает переполнение 16-битных операндов. Пришлось ограничится простым вариантом с 10-ю сложениями операндов, и 8-ю сложениями результата. В этом варианте операнды удваиваются и максимальная глубина без переполнения – 5. Правда в конце оптимальная по скорости глубина оказалась 4.
5. В конце я обнаружил, что GCC 4.0 плохо оптимизирует свои intrinsic SSE функции – оставляет некоторое количество бессмысленных копирований регистров. Переписал на ассемблере единственную функцию умножения 4-х строк и получил последнее ускорение.

Решение тов. alexBlack

Для эффективного использования кэша процессора матрица разбивается на матрицы меньшего размера. В [1] приводится вариант разбиения с изменением индекса. У меня использован другой вариант. Каждая меньшая матрица физически расположена в непрерывном участке памяти и все три матрицы представлены как матрицы, состоящие из M2xM2 матриц меньшего размера M1xM1:

type
tm1 = array [0..M1-1, 0..M1-1] of smallint;
tm2 = array [0..M2-1, 0..M2-1] of tm1;

Это объявление типа для исходных матриц A и B. Для матрицы результата объявление такое-же, но тип элементов integer. Естественно, что теперь размер исходных матриц может принимать только ряд дискретных значений и с лишними частями нужно что-то делать. У меня они просто заполняются нулями и используются в вычислениях.

Что касается констант, то после нескольких экспериментов было выбрано M1 = 64 и M2 = 79 (Все три меньших матрицы помещаются в кэш первого уровня, а 64*79=5056, что покрывает максимальный размер матрицы). Такое разбиение не очень сильно усложняет ввод/вывод матриц.

Теперь о вычислениях. Математически формула умножения матриц ничем не отличается для элементов, которые сами являются матрицами. Отсюда код:

procedure MulM1xM1(var R:tmI1; var A, B:tm1);
var i, j, k:integer;
s0 : integer;
begin
for i := 0 to M1-1 do begin
for j := 0 to M1-1 do begin
s0 := 0;
k := 0;
while k < M1 do begin
s0 := s0 + A[i,k+ 0]*B[j,k+ 0]+
// .... ....
A[i,k+15]*B[j,k+15];
inc(k, 16)
end;
R[i,j] := R[i,j] + s0;
end;
end;
end;

procedure MulMatrix;
var i, j, k:integer;
begin
for i := 0 to M-1 do begin
for j := 0 to M-1 do begin
for k := 0 to M-1 do begin
MulM1xM1(R[i, j], A[i, k], B[j, k]);
end;
end;
end;
end;

Здесь M=[N/M1] - количество маленьких матриц в строке большой матрицы для заданного размера матрицы N.

Обратите внимание, для матрицы B изменен порядок следования индексов и, соответственно, используется транспонированная матрица. (Транспонирование производится при чтении и не сказывается на общем времени вычислений). Такой порядок индексов позволяет читать элементы массива по строкам, что значительно эффективнее. Кроме того, самый внутренний цикл частично развернут. Это еще один из способов оптимизации за счет устранения точек ветвления.

Замечу, что прямая реализация на ассемблере ничего не дает в плане ускорения работы. Delphi генерирует достаточно эффективный код. В частности, при обращении к элементам массива не используется умножение, а делается что-то что-то вроде:

movsx eax, [esi+ebx*2+0]
movsx edx, [edi+ebx*2+0]
imul eax, edx
;...

здесь ebx - счетчик k внутреннего цикла, 2-ка - размер элемента (элементы матрицы хранятся как signed word).

Далее самый внутренний цикл (while k) заменен ассемблерной вставкой с использованием команд MMX. На самом деле была заменена вся процедура, но здесь я приведу только код внутреннего цикла:

asm
mov esi, PA // @A[i, 0]
mov edi, PB // @B[j, 0]

movdqa xmm0, dqword ptr [esi+ 00]
movdqa xmm1, dqword ptr [esi+ 16]
// и т.д. для xmm2-xmm6 .... для 32, 48, 64, 80, 96
movdqa xmm7, dqword ptr [esi+112]

PMADDWD xmm0, dqword ptr [edi+ 00]
PMADDWD xmm1, dqword ptr [edi+ 16]
// и т.д. для xmm2-xmm6 .... для 32, 48, 64, 80, 96
PMADDWD xmm7, dqword ptr [edi+112]

paddd xmm0, xmm1
paddd xmm0, xmm2
// .... для xmm3 .. xmm7
paddd xmm0, xmm7

pshufd xmm2, xmm0,1110b
paddd xmm0, xmm2
pshufd xmm2, xmm0,0001b
paddd xmm0, xmm2

movd s0, xmm0
end;

Как видим этот блок вычисляет произведение одной строки матрицы 64x64 на строку другой матрицы 64x64.

При реализации столкнулся с тем, что в Delphi нет директивы выравнивания по границе параграфа. Максимум Align8. Пришлось выравнивать вручную добавлением лишних переменных.

Последние команды - сложение по горизонтали можно заменить командами PHADDD:

// Сложение по горизонтали
pxor xmm1, xmm1
//PHADDD xmm0, xmm1
db $66, $0F, $38, $02, $C1
//PHADDD xmm0, xmm1
db $66, $0F, $38, $02, $C1

К сожалению, Delphi о них не знает, да и время выполнения почти не меняется. Кстати, эти команды отказался выполнять имеющийся у меня (не самый старый) процессор AMD.

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

asm
PADDD xmm0, dqword ptr [ebx]
movdqa dqword ptr [ebx], xmm0
end

Не буду приводить полный код процедуры, т.к. он достаточно большой, а основные части его уже приводились.

Архив конкурса
(цитата)

Все желающие могут скачать архив конкурса. В нём содержится:

• простая тестирующая система, взятая с олимпиад по программированию;
• генератор тестов и правильных ответов к ним;
• программы всех участников (только их последние версии);
• исходные коды тех участников, которые пожелали ими поделиться (мой код прилагается);

К сожалению, победитель не даст исходный код, но он обещал позже поделиться основными идеями. Когда он это сделает, я напишу здесь же или в комментариях. [ Дополнение 23.02.2010: тов. venco описал своё решение в комментариях, по-моему, всё понятно, нужно попытаться повторить результат ].

Скачать архив (2 Мб).

Поясню, как этим пользоваться. Программа tchoose написана не мной, о истории её создания можно узнать в файле readme.txt. Программа вычисляет процессорное время, но с погрешностью где-то в полсекунды. Когда мне казалось, что какие-то программы участников конкурса работали на одном из тестов слишком долго (дольше, чем ожидалось, исходя из уже пройденных тестов) я всегда перетестировал данный тест трижды и выбирал лучшее время. Но поскольку почти у всех участников программы сильно отличаются по времени и проблем с лидерством не возникло, то эта погрешность не повлияла на результаты конкурса. Система, кроме запуска Вашей программы, также проверяет ответ на правильность, но это происходит гораздо дольше, чем работа самой программы : ) (такой чекер), поэтому реальная проверка одного решения занимает много времени.

Для начала нужно зайти в папку tchoose/tests и сгенерировать тесты программой gen.exe (исходник этой программы там же). После этого нужно запустить solve.bat — этот скрипт решает все тесты с помощью моей программы solution.exe (это моя лучшая кубическая версия, 14,5 с., правильность гарантирована), и сохраняет ответы. После этого в директории появятся файлы 01, 01.a, 02, 02.a и т. д. (файлы 00 и 00.a не генерируются, они есть там с самого начала). Будте внимательны, все тесты с решениями занимают 4400 Мб, убедитесь, что на вашем диске хватает места.

После того, как тесты сгенерированы, зайдите в папку tchoose/solutions и положите туда свою программу (программы всех авторов с исходниками хранятся там же). Теперь запустите tchoose.exe и следуйте интуитивно понятному интерфейсу (нажмите внизу кнопку проверки решений). Если возникнут проблемы, нужно прочитать readme.txt, в котором также написано, как можно изменить файлы, ограничения по времени и по памяти.

Можете запустить все программы на всех тестах и убедиться в том, что окончательная таблица рейтинга составлена правильно. Конечно, времена работы у Вас могут отличаться.

В директории с программой есть файл check.cpp — это чекер, написанный по правилам этой системы. Работает на основе кода из testlib.h (этот файл там же). Чекер простой: сверяет два файла, игнорируя различия в пробелах, табуляции и пр.

Если по какой-то причине tchoose у Вас не работает, то я ничего не могу сделать — у Вас проблемы с Windows. Любителей Linux прошу меня простить, аналогичной программы для Linux я не знаю.