Перемножение матриц в С++

[arseniiv]

Давайте всё-таки использовать [syntax lang="cpp"]…[/syntax] или там [syntax lang="c"]…[/syntax], раз уж они есть.

[Andrey_Kireew]

хорошо!

[mihaild]

(Оффтоп)

Dmitriy40 в сообщении #1225404 писал(а):

скорость работы с памятью считаем уже оптимизированной и не тормозящей процесс

На практике перемножение матриц в несколько потоков упирается в чтение из памяти обычно.

Andrey_Kireew в сообщении #1225408 писал(а):

кстати экспонента алгоритма перемножения matlab примерно 2.69, это лучше, чем алгоритм Штрассена

А как вы это меряете? Чтобы поймать разницу в $0.04$ надо в показателе надо мерять очень аккуратно (и на больших матрицах, чтобы перестало быть важным, какая доля влезает в кеш какого уровня)

Andrey_Kireew в сообщении #1225475 писал(а):

теперь цикл выполняется за 1.5 секунды

А какое железо? (главное - какие SIMD инструкции поддерживает процессор)
Вы же с оптимизацией собираете? (просто полторы секунды - это всё равно очень много; либо железо старое, либо что-то не так)

[Andrey_Kireew]

Я посчитал время для матрицы 10000x10000 получилось 46.5 секунд кажется, ну и из пропорции
$46.5\cdot 1000^x=0.1$\cdot 10000^x
$x=log_{10}{465}$ (это по памяти, может быть не совсем точно, те расчёты я не сохранил)

я проверяю на ноте i5-3317u (2 ядра, 2.4 ГГц в режиме турбо), но не суть важно, есть железо и по мощнее, просто это самый трудоёмкий фрагмент всей моей программы, причём он далеко не первостепенной важности, т.е. не вписывается он в приложение так как есть, его или оптимизировать, или убирать вообще

но 1.5 секунды это уже терпимо, попробую ещё оптимизировать по Штрассену и уменьшить общее количество циклов

[Xaositect]

Andrey_Kireew в сообщении #1225503 писал(а):

Я посчитал время для матрицы 10000x10000 получилось 46.5 секунд кажется, ну и из пропорции
$46.5\cdot 1000^x=0.1$\cdot 10000^x
$x=log_{10}{465}$ (это по памяти, может быть не совсем точно, те расчёты я не сохранил)

Такие вещи просто так нельзя считать по двум точкам. Выбираете достаточно много достаточно больших значений $n$, для каждого $n$ делаете несколько замеров умножения матриц $n\times n$, замеряете время $t$ и строите регрессию $t = C n^{\omega}$.

[mihaild]

Andrey_Kireew в сообщении #1225503 писал(а):

но 1.5 секунды это уже терпимо, попробую ещё оптимизировать по Штрассену и уменьшить общее количество циклов

Попробуйте лучше взять готовую библиотеку (если есть возможность - MKL). Пытаться обогнать готовые решения в таких задачах - почти всегда пустая трата времени.

[Andrey_Kireew]

:), да я обогнать и не пытаюсь, но если есть возможность - стараюсь обходится без сторонних библиотек, чтобы в них не запутаться, да и не надо вникать в интерфейсы новых непонятных функций и под них подстраиваться ...

Честно говоря до вчерашнего дня я в перемножении матриц вообще не видел никакой проблемы.

Надо подумать, может быть и библиотеку подключу, только MKL наверное не получится, там онлайн - инсталятор и всё такое... Может посоветуете что по проще?

[mihaild]

openblass
Тот же API. Работает помедленнее, но не очень сильно.

[YYSS]

Andrey_Kireew в сообщении #1225268 писал(а):

В моё приложении на С++ основное время затрачивается на перемножение квадратных матриц...

В этом случае лучше использовать специализированные библиотеки линейной алгебры: у меня на 4-ех ядерном процессоре, поддерживающем команды fma, например, матрицы размерами 10000х10000 библиотека Intel MKL перемножает за 10 секунд. Т.е. матрицы размерами 1000х1000 будут перемножаться примерно за 0.01 с. или в 1000 раз быстрей, чем у вас.

[bondkim137]

Во-первых, читать столбцы - не очень хорошо. Из-за особенностей работы кешей. Следует за квадратичную сложность сначала развернуть одну матрицу, причем так: читаем строку, пишем столбец.

Потом перемножать строка на строку. Внутренний цикл перевернуть и выкинуть оттуда все, кроме адресации по индексу и одновременно счетчику, вроде такого:

Используется синтаксис C++

for (int j = N - 1; j >= 0; j--)
  s+= pA[j]*pB[j];
 

И чего точно не стоит делать, это сохранять промежуточную сумму в память.

Затем, если написать этот цикл на ассемблере с использованием SIMD и, при необходимости, в несколько потоков, то все вопросы окончательно отпадут.

[YYSS]

bondkim137 в сообщении #1225733 писал(а):

Потом перемножать строка на строку...

Блоки строк.

[Dmitriy40]

Можно обойтись без перевора всей матрицы, пересохранять лишь текущий столбец в дополнительный массив. Сложность та же, дополнительной памяти в $n$ раз меньше, исходная матрица не изменяется. Параллелить мне кажется логичнее самый внешний цикл (меньше накладных расходов).
Проше и логичнее воспользоваться готовой библиотекой, как выше указано, она достигает практически пиковой скорости вычислений.

[bondkim137]

Dmitriy40 в сообщении #1225755 писал(а):

пересохранять лишь текущий столбец в дополнительный массив

Плюсую, как я сам не догадался :)

Dmitriy40 в сообщении #1225755 писал(а):

Параллелить мне кажется логичнее самый внешний цикл

Его и надо. Если вообще надо.

Dmitriy40 в сообщении #1225755 писал(а):

Проше и логичнее воспользоваться готовой библиотекой

Это если, конечно, спортивный/учебный интерес к низкоуровневой оптимизации отсутвует

[Dmitriy40]

Если кому-то захочется повозиться с оптимизацией именно на языке С, то вот неплохой пример как это уже проделывали. Что характерно, после всех усилий (причём пошедших дальше чем здесь предлагали выше) библиотеку Intel MKL даже близко не догнали по скорости, разница осталась примерно порядка (10-15 раз)! Что неплохо иллюстрирует разницу между практикой и учёбой, если нужен результат на практике, то намного выгоднее воспользоваться плодом труда коллективного разума в виде библиотеки.

-- 15.06.2017, 22:11 --

Добавлю ещё две полезные в учебных целях ссылки на другой форум про оптимизацию перемножения матриц: раз (на C), два (ахтунг, ассемблер!).

[YYSS]

Dmitriy40 в сообщении #1225844 писал(а):

если нужен результат на практике, то намного выгоднее воспользоваться плодом труда коллективного разума в виде библиотеки.

Не всегда: например, для перемножения матриц
https://software.intel.com/en-us/forums ... pic/288850
а это для диагонализации
https://software.intel.com/en-us/forums ... pic/288316
https://software.intel.com/en-us/forums ... pic/288218
https://software.intel.com/en-us/forums ... pic/287728
https://software.intel.com/en-us/forums ... pic/289711
https://software.intel.com/en-us/forums ... pic/290238