Итак, конкурс по перемножению матриц, в котором нужно было перемножить две матрицы размером 5000×5000, завершён. Спасибо всем участникам, благодаря им я узнал много нового: конкурс заставил меня делать такие вещи, в которые я сначала не верил или в которых не разбирался. Здесь я напишу несколько итоговых комментариев и выложу архив с решениями и тестами.

Окончательный рейтинг

В конкурсе приняли участие 8 человек. С большим отрывом от остальных победил тов. venco, второе место занимают сразу два участника: blackmirror и я (погрешность измерения времени оказалась чуть больше, чем разница между работой наших программ на самых больших тестах).

  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 с.

Архив конкурса

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

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

К сожалению, победитель не даст исходный код, но он обещал позже поделиться основными идеями. Когда он это сделает, я напишу здесь же или в комментариях. [ Дополнение 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 я не знаю.

Моё решение (кратко)

Во-первых, сначала я написал самую тупую версию, которая умножает строку на столбец в тройном цикле (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, и это совершенно ничего не изменило. Даже чуть медленнее стало. Но до конца конкурса оставалось около часа, и я так и не успел разобраться до конца.

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

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

Пояснения к конкурсу

Пояснения к правилам

Поясню, какова была моя цель и причина указанных ограничений. Размеры матрицы, конечно же, были выбраны наугад, на самом деле можно было выбрать и другие, но 4000 — маловато, а 6000 — «неровно», поэтому 5000 взято как просто «удобное» число. Умножая время работы на 8 (или на 7, в зависимости от алгоритма), можно приблизительно оценить работу программы на порядках матрицы 10000 и 20000 (при условии наличия нужного объема памяти). Почему выбран 1Gib? Чтобы реализация методов типа Штрассена совсем «в лоб» не проходила. Откуда взялось 600 в ограничениях на входные данные? Я хотел, чтобы ответ помещался в 32 бита со знаком. Конечно, эти ограничения можно было сделать чуть другими, скажем, были бы входные данные — двойные слова, а на выходе — четверные. Не суть важно, конечно pmaddwd уже не работала бы, но нашлась бы другая удобная инструкция.

Самый главный вопрос: почему нельзя распараллелить? И правда, многие делали разумные замечания: «используй CUDA и т. п., будет быстрее». Дело в том, что перемножение матриц в моих задачах является частью уже распараллеленной программы, каждый поток которой должен много раз перемножать матрицы, поэтому распараллеливать нечего. Далее, если бы нужно было решать перемножение матриц как отдельную задачу, то я бы решал её на кластере за гораздо меньшее время. Почему решено делать ввод данных из файла? Многие из тех задач, которые я решаю, должны считывать данные из файла, но не 250 Мб, а гораздо больше. Поэтому было решено, что в этом конкурсе будет соревнование не только на вычисление, но и на ввод / вывод данных.

Смысл конкурса

Вы, конечно же, знаете, что задача перемножения матриц изучается очень давно, и для её решения существует множество разнообразных алгоритмов. Многие «теоретики» в своих научных статьях часто упоминают сложность работы каких-то алгоритмов, редко задумываясь об их практической полезности. В рамках конкурса я хотел понять, на что именно можно рассчитывать на практике. До определённого момента я не верил, что программа может работать сильно быстрее минуты (две недели назад я не знал границ применимости SSE, а компилятор Intel с опцией /QxSSE4.1 давал как раз программу, работающую одну минуту (без оптимизации по блокам)). С другой стороны, я знаю, что многие любят пользоваться стандартными библиотеками (типа BLAS или MKL), получая при этом не слишком быстрые программы, но все-таки достаточно оптимальные. Все такие библиотеки, расчитанные на общий случай, по идее не должны даже рядом лежать с ручной реализацией, которая разрабатывалась для частного случая. Это мне тоже нужно было показать. Но тов. ghoul меня сильно удивил, когда показал программу, написанную с использованием BLAS, которая работала гораздо быстрее моих первых версий. С этого момента мне стало понятно, что из 30 секунд можно выйти (потом я проверял MKL, она выдавала около 25 секунд). Если бы я работал над этой задачей один, вне конкурса, то никогда бы не додумался, как можно уложиться в 10 секунд. А когда участники постоянно присылали программы, работающие все быстрее, мне ничего не оставалось, как не отставать от них.

В результате конкурса было показано, что при таких ограничениях матрицы можно перемножать гораздо быстрее, чем это умеют делать функции стандартных библиотек. Также было показано, какого именно результата можно ожидать на практике. Экстраполируя результаты, можно делать выводы. Например, если входные данные будут больше (скажем, 32 бита и целые), то ожидаемое время работы составит 14-16 с. В любом случае, полученные данные будут полезны не только тем, кому нужно решать такую задачу быстро, но и тем, кто, занимаясь голой теорией, забывает о реальной практической значимости своих результатов. Например, многие (причём уважаемые люди) мне постоянно утверждают, что уже давно придуман алгоритм, работающий со сложностью O(n2.36) и что это «самый лучший» результат. Теперь я могу сказать таким людям, что если они могут написать этот алгоритм так, чтобы оказаться в таблице с рейтингом на первом месте, то замечательно, а пока это только слова, не имеющие никакой ценности при решении практических задач.

Личное отношение к задаче

Здесь я скажу пару слов о самой задаче и выражу своё личное мнение.

Итак, некоторые люди на форумах говорили, что я не прав, назвав задачу простой. Дело в том, что я делю задачи на простые и трудные (или труднорешаемые). Почти все полиномиально-разрешимые задачи являются простыми. В практических приложениях сложность простых задач легко «придавить» сверху количеством ядер. Например, то же перемножение матриц можно решать при размерах до 10000 всего на 8 ядрах за то же время. То есть можно удвоить размерность, при этом число ядер тоже увеличится в определённое число раз. В труднорешаемых задачах с экспоненциальной сложностью умножение числа ядер на константу нужно делать при увеличении размера задачи всего на единицу. А про удвоение размера задачи вообще говорить страшно. Таким образом, задавить сложность числом узлов в кластере не получится и выкручиваться нужно уже не программным, а больше алгоритмическим способом. Вот мы ускорили перемножение матриц более чем в 100 раз, а в труднорешаемых задачах (на репрезентативных размерностях) ускорять приходится минимум в миллионы раз, и это только на начальной стадии. Как Вы думаете, что все-таки сложнее?

Таким образом, сама по себе задача перемножения матриц является простой. Это не уменьшает заслуг тех, кто придумал для нее замечательные алгоритмы и тех, кто умеет хорошо программировать, но в настоящее время имеются гораздо более важные задачи, почти 100% которых до сих пор не решены не только эффективно, но даже теоретически для них не получен результат. Пример? Ну возмите самую известную задачу о расстановке n ферзей на шахматной доске nxn. В прошлом году посчитали ответ для n=26. Сейчас по этой задаче даже неизвестна асимптотика роста числа расстановок с ростом числа n (хотя есть гипотеза, что это примерно n!/2.54n). А улучшать простую задачу можно сколько угодно, даже для чудовищных размеров нужно (и почти всегда можно) просто подобрать достаточное количество ядер. В любом случае, я считаю, что решать трудные задачи намного интереснее, хотя сам когда-то занимался только полиномиальными (например, мне очень нравилась задача о максимальном потоке, транспортная задача и прочие задачи оптимизации, разрешимые в классе P).

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

Всем огромное спасибо за помощь!

Обсуждение конкурса на форуме в этой теме.