Конкурс по обращению матрицы завершён. Победителем объявляется неоднократный участник моих конкурсов alexBlack.

  1. alexBlack :: Delphi :: 267,8 с.
  2. зщл :: g++ 4.6.2 :: 836,0 с.
  3. Sher :: Visual C++ / Asm :: 1001,5 с.
  4. Student :: Visual C++ :: 1691,8 с.

Поздравляем победителя. Должен отметить, что отрыв от остальных участников достаточно солидный для такой задачи.

Думаю, давать систему тестов нет смысла, генерировались они всё равно случайно.

Моё решение

Свою более-менее быструю программу я решил написать только после того, как победитель прислал свою версию, мне стало очень интересно, как получилась такая высокая скорость. Моя прямолинейная реализация работала чуть меньше 2400 секунд, она была написана только чтобы подготовить тесты. Последняя версия, но не самая быстрая, работала 231 секунду. Опишу кратко, как добиваться такого ускорения.

  • Пишется блочный метод Гаусса, когда матрица разбивается на блоки, умещающиеся в кэш L2. Даже «на коленках» написанная блочная версия уже даёт 2095 секунд. В ней обращение матрицы размером 2n × 2n заменяется на 5 умножений и 2 обращения матриц размером n×n. То есть сложность алгоритма уже чуть меньше, чем кубическая. Подробное описание этого алгоритма есть в книге Гантмахер Ф. Р. «Теория матриц» (книга свободно валяется в сети) в главе II, §5 «Разбиение матрицы на блоки. Техника оперирования с блочными матрицами. Обобщенный алгоритм Гаусса». Подробнее этот алгоритм я опишу в другой раз.
  • Профилировщик сказал, что 99% времени работы программы занимает как раз перемножение матриц. Переписываем его на ассемблере (пока без SSE, используя только div для взятия остатка), аккуратно распределяя регистры. Получается 1165 секунд.
  • Следующий шаг — нужно убрать взятие остатка в самом глубоком цикле перемножения матриц. Если мы считаем произведение C=A⋅B, то число Cij получается как сумма произведений Aik ⋅ Bjk (т. к. матрица B, естественно, транспонирована), так вот, эту сумму можно сразу брать по модулю, а можно потом, после завершения цикла по k. Промежуточный результат укладывается в 96 бит, а затем остаток вычисляется двумя операциями div. Таким образом, вместо n3 делений нужно всего 2n2. Это дало мне 415 секунд.
  • В методе Гаусса деление на 231-1 все равно остаётся. Многие знают, что делить на степень двойки легко и приятно, однако не многие в курсе, что деление на степень двойки без единицы тоже можно выполнить битовыми трюками. Подробнее об этом я напишу в ближайших постах. К сожалению, эта оптимизация ничего не дала, так как блочный Гаусс и трюк с int96 уже ускоряли все достаточно сильно. Метод Гаусса на маленьких блоках все равно работал быстро что до применения трюков с битами, что после.
  • Переделываем перемножение матриц на SSE. Первая попытка замедлила код в 2 раза, но потом я внимательнее перечитал sse instruction set, и понял, что неаккуратно подобрал команды. Правильная подборка и последовательность команд movdqa, shufps, pmuldq и pextrd дала 231 секунду. Здесь также использовался тот факт, что 4 числа, ограниченных 231-1 в сумме произведений еще влезают в int64, добавлять сумму в int96 нужно только каждый четвертый раз.
  • Я не пробовал учесть особенностей кэша L1 при перемножении матриц. Уверен, что с ним можно довести до 200 секунд. Кроме того, программа моя написана не очень аккуратно (поэтому я её не покажу), так что полученное время — это ещё не предел.

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

Всем спасибо за участие! Приходите ещё.

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