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

[Lia]

YYSS заблокирован как злостный клон.

[Andrey_Kireew]

dsge в сообщении #1225391 писал(а):

Вы можете скомпилировать м-файл и потом посмотреть С++ коды, какой получится алгоритм

Всё таки решил последовать Вашему совету и посмотреть,
вот что получается после matlab-coder

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

14   /* Function Definitions */
   15   void PROD(const real_T A[10000], const real_T B[10000], real_T C[10000])
   16   {
   17     int32_T i0;
   18     int32_T i1;
   19     int32_T i2;
   20     for (i0 = 0; i0 < 100; i0++) {
   21       for (i1 = 0; i1 < 100; i1++) {
   22         C[i0 + 100 * i1] = 0.0;
   23         for (i2 = 0; i2 < 100; i2++) {
   24           C[i0 + 100 * i1] += A[i0 + 100 * i2] * B[i2 + 100 * i1];
   25         }
   26       }
   27     }
   28   }
   29   
   30   /* End of code generation (PROD.c) */
 

в этом коде даже порядок циклов не оптимален,
если же создавать не .DLL а .mex функцию, то работает быстро.
Cтал разбираться - во втором случае С - код уже другой, и его намного больше,
для перемножения там используется dgemm() из библиотеки BLAS.h,
сама эта библиотека оказывается есть среди файлов matlab.

Попробовал использовать dgemm() в своей программе, как ни странно всё получилось,
теперь перемножение выполняется даже немного быстрее чем в matlab

[Andrey_Kireew]

В дополнение к теме. Matlab поддерживает CUDA вычисления. Я попробовал, на моей простенькой видеокарте NVIDIA GT1030 - перемножение матриц выполняется в 50 раз выстрее чем на ЦПУ (intel core i7 6700). Матрицы формата double 10000х10000 (больше не понещается в видеопамять, но думаю это не является большой преградой, т.к. матрицы можно разбивать на блоки).

[photon]

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

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

   for (k=kk; k<kk+nb; k++) 
       val+=*(BT+k+i*NStr)**(M+k+j*NStr);

 

Это можно еще улучшить (не проверял, но не уверен, что компилятор всегда соптимизирует такие конструкции):
1) j*NStr можно вычислить за циклом, поскольку оно в нем не меняется
2) переменную k можно вообще выкинуть

сделать что-то типа:

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

auto it1 = BT + kk + i * NStr;
auto it2 = M + kk + j * NStr;
auto end = it1 + nb;
for (; it1 < end;)
{
    val += *it1++ * *it2++;
}
 

[bondkim137]

- Сравнение с числом (end) не бесплатно. Сравнение с нулем после уменьшения счетчика бесплатно. Потому циклы обычно разворачивают от N до нуля.
- Инкрементировать указатель, как правило, не бесплатно. Разыменовывать по индексу, как правило, бесплатно.
если не расчитывать, что компилятор сильно умный, то лучше так:

for (int i = (end - it1) - 1; i >= 0; i--)
{
    val += it1[i]*it2[i];
}

[Dmitriy40]

Сами же сказали что бесплатно сравнение с нулём после уменьшения, в коде же сначала сравнение, потом уменьшение. Лучше уж так:

for (int i = (end - it1); --i >= 0; ) { ... }

Кроме того, оба этих утверждения

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

Сравнение с числом (end) не бесплатно. Сравнение с нулем после уменьшения счетчика бесплатно.

недостаточно обоснованы (мягко говоря), я специально проверил вот такой код инструментом Intel® Architecture Code Analyzer

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

@1:     dec     ECX
        jge     @1
@2:     cmp     ECX,100
        jl      @2

на архитектуре Nehalem первая команда перехода вообще не спаривается, на следующих архитектурах (Sandy Bridge, Ivy Bridge, Haswell, Broadwell, Skylake) оба условных перехода спариваются совершенно одинаково. От разрядности (x32 или x64) не зависит. От типа числа для сравнения (регистр, непосредственный операнд, прямо или косвенно адресуемая ячейка памяти) тоже не зависит. Т.е. оба процитированных утверждения неверны, по крайней мере для перечисленных архитектур Intel.
Потому сравнение с конечным числом выгоднее счётчика на старых процессорах и без разницы на относительно новых. И соответственно проходить цикл от N к 0 или от 0 к N для распространённых процессоров разницы нет (если забыть про необходимость дополнительного inc ECX во втором случае или сравнивать не счётчик цикла, а указатели, которые и так надо модифицировать).

[bondkim137]

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

Сами же сказали что бесплатно сравнение с нулём после уменьшения, в коде же сначала сравнение, потом уменьшение. Лучше уж так

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

-- 01.12.2017, 15:14 --

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

недостаточно обоснованы (мягко говоря), я специально проверил вот такой код инструментом Intel® Architecture Code Analyzer

У вас там, если я правильно понял, бесконечный цикл во втором случае. А если вы инкремент указателя/счетчика вставите и получите три команды, оно может как спариться (обойтись бесплатно), так и не спариться - в т.ч. зависит от железа. Также может спариться, но занять место в конвеере, где что-нить другое могло бы спариться (вы же не пустые циклы гонять собираетесь) - а в первом случае декремент скорее всего спарится с чем-нить полезным.

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

[Dmitriy40]

bondkim137
Я же добавил про инкремент счётчика, от него можно отказаться используя указатель в качестве счётчика циклов (сравнивая указатель с конечным его значением). В этом случае всё будет бесплатно. Да ещё и регистр счётчика цикла сэкономится.
Насчёт "своего железа", я же написал на каких архитектурах это верно, а не только на моём железе.

Добавлю, оба варианта спариваются и не спариваются совершенно одинаково (начиная с Sandy Bridge и более поздних), потому слова о возможном неспаривании вообще роли не играют, оба варианта при этом будут (не)спарены совершенно одинаково.
Потому в чём выгодность именно уменьшения счётчика до нуля всё ещё не очень понятно. Экономится лишь команда inc, от которой можно вообще отказаться вместе со счётчиком цикла и сэкономить регистр.
Да и вообще это уже поиск блох, имхо, оптимизировать надо в другом месте.

[bondkim137]

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

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

Пардон, но этот указатель нужно же еще инкрементить.
Еще раз:

Вариант прямой:
1) разыменовываем указатель
2) делаем полезную работу
3) инкрементим указатель
4) сравниваем указатель с ограничителем
5) условный переход

Вариант обратный:
1) разыменовывем указатель+индекс (очень давно это везде одна команда и индексирование является бесплатным)
2) делаем полезную работу
3) декрементим индекс
4) делаем условный переход

Совершенно из общих соображений обратный вариант на одну команду короче.
А если вспомнить, что есть команды вроде loop, которые декрементят счетчик и выполняют условный переход сразу, то вообще так:

1) разыменовывем указатель+индекс (очень давно это везде одна команда и индексирование является бесплатным)
2) делаем полезную работу
3) декрементим индекс+делаем условный переход

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

(Оффтоп)

loop на практике как правило, не нужен, т.к. с большой вероятностью декремент указателя удастся запараллелить с полезной работой. Особенно, если последней командой блока полезной работы будет не меняющая флаги команда (в т.ч. весь SIMD) - тогда на всякий случай рекомендовано поднять декремент на одну/несколько позиций выше (компиляторы это тоже делают как правило). Так мы делали с времен появления первого Pentium и до недавних времен, пока я еще занимался низкоуровневыми оптимизациями.

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

[Dmitriy40]

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

А если вспомнить, что есть команды вроде loop, которые декрементят счетчик и выполняют условный переход сразу,

Угу, команда удобная, вот только сильно медленная (и да, я её только что протестировал). За подробностями к Агнеру Фогу или тестам.

Итого.
Я согласен с Вами что прямой ход (с 0 до N) на одну команду длиннее и на такт дольше (не всегда, в зависимости от условий).
Я не согласен что

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

Сравнение с числом (end) не бесплатно. Сравнение с нулем после уменьшения счетчика бесплатно.

Как я привёл данные выше, именно в такой формулировке это не так, и то и другое давно одинаково бесплатно (для современных процессоров Intel). Если не брать остальную часть цикла, а только процитированное! Для всего цикла в целом картина вполне может измениться. Т.е. эта формулировка недостаточно корректна.

[bondkim137]

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

Угу, команда удобная, вот только сильно медленная (и да, я её только что протестировал). За подробностями к Агнеру Фогу или тестам

Зачем к Фогу сразу? За мою карьеру с 386-х она периодически то нормально работает, то медленно :) Пару раз быстро (быстрее двух команд) работала. Потому для универсальности ее обычно не рекомендуют использовать, оговаривался выше - если грубо вслепую писать, то лучше просто декремент засунуть повыше в тело цикла, тогда он почти наверняка почти везде параллельно с чем-нить выполнится. Так мы всегда и делали.

Любопытно: на сколько тактов сейчас у вас медленнее loop получился (и на каком процессоре), чем декремент + условный переход?

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

Как я привёл данные выше, именно в такой формулировке это не так

Вы правы. Впрочем, я все-таки имел в виду определенный конкретный контекст :)

[Dmitriy40]

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

Любопытно: на сколько тактов сейчас у вас медленнее получилось (и на каком процессоре) чем декремент + условный переход?

Реальные тесты я провожу на Haswell (за другие компы идти лень, да и там тормозные ноуты), декремент с условным переходом примерно 1 такт занимает, loop занимает примерно 4 такта. Это с пустым циклом. Небольшие отклонения (процентов 5) связаны с погрешностью измерений.

[bondkim137]

(Оффтоп)

понятно, спасибо

[GAA]

Dmitriy40, почему dec/inc?

Intel® 64 and IA-32 Architectures Optimization Reference Manual, 2016 (а также в более ранних редакциях) писал(а):

Assembly/Compiler Coding Rule 33. (M impact, H generality) INC and DEC instructions should be replaced with ADD or SUB instructions, because ADD and SUB overwrite all flags, whereas INC and DEC do not, therefore creating false dependencies on earlier instructions that set the flags.

В реальном коде это сопли и тестировать очень трудно, а рассчитать не просто. Но почему не прислушаться к рекомендациям Intel? И без всяких тестов использовать sub/add? (Это я о 64b. 32b не тестил уже давно.)

-- Сб 02.12.2017 13:29:42 --

Аналогично test vs cmp

Intel® 64 and IA-32 Architectures Optimization Reference Manual писал(а):

Assembly/Compiler Coding Rule 19. (M impact, ML generality) Employ macro-fusion where possible using instruction pairs that support macro-fusion. Prefer TEST over CMP if possible. Use unsigned variables and unsigned jumps when possible. Try to logically verify that a variable is nonnegative at the time of comparison. Avoid CMP or TEST of MEM-IMM flavor when possible. However, do not add other instructions to avoid using the MEM-IMM flavor.

[Dmitriy40]

GAA в сообщении #1271021 писал(а):

Но почему не прислушаться к рекомендациям Intel? И без всяких тестов использовать sub/add?

Потому что в рамках обсуждаемого вопроса (спаривание условного перехода и предыдущей команды сравнения) это не имеет значения, они спариваются одинаково. И занимают одинаковое время.
Кроме того, в процитированных словах прямо сказано что возникают ложные зависимости к близкой следующей команде модифицирующей (устанавливающей) флаги, команда условного перехода к таким не относится, потому эта оптимизация на данную связку команд не влияет.

Насчёт TEST vs CMP, да, это знаю, но CMP значительно удобнее человеку т.к. не требует для завершения точного равенства, достаточно проверить на (не) превышение порога. В ситуациях когда начальное смещение неизвестно или приращение не является делителем общего размера (а иногда это заранее неизвестно) TEST неприменима. Ну и чтоб не думать, удобно везде ставить CMP вместо TEST для сравнения чисел, а TEST оставить для двоичных полей и логических типов. Ну а конкретно здесь и CMP прекрасно спаривается, так что замена на TEST смысла не несёт. Привычка же заменять всегда не вырабатывается из-за неполной взаимозаменяемости TEST и CMP (далеко неполной, вот недавно писал цикл, так самым удобным было сравнить два указателя не на равенство, а на "переворачивание", т.е. когда меньший становится точно больше большего, причём точного совпадения их могло вообще не быть, для чётного количества элементов). Ну и обычно на такие мелочи просто не обращаю внимания, это ж не основная нагрузка, а какие-то жалкие единицы процентов скорости.

Ещё раз подчеркну, с моей стороны здесь речь шла об очень узком смысле, исключительно о спаривании условного перехода и предыдущего сравнения. В более общем/универсальном случае рекомендации Intel безусловно верные.

(Исключение)

Хотя есть и исключение, команда inc/dec занимает на байт меньше и если затык именно в скорости выборки и декодирования инструкций (а для AVX кода это нередкий случай, команды длинные) может быть быстрее чем add/sub, на длинных циклах, не помещающихся в кэши декодированных инструкций. Правда тогда эти единицы тактов роли уже и не играют, но формально таки быстрее ...

-- 02.12.2017, 16:16 --

(Пара примеров TEST vs CMP, для иллюстрации)

Из недавнего:

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

;Суммирование массива попарно с начала и конца
@sum1:  fld     extended [EDX]
        fld     extended [ECX]
        faddp
        fstp    extended [EDX]
        add     EDX,10
        sub     ECX,10
        cmp     EDX,ECX         ;TEST неприменима т.к. совпадения указателей может и не быть, а корректно вычислять количество элементов и для чётного и для нечётного случая поленился, да и регистр для счётчика жалко
        jc      @sum1

;Суммирование элементов с пропусками
        mov     ECX,10*101      ;101 элемент в массиве
@sum2:  fadd    extended [ESI]
        add     ESI,10*1        ;Поставив любой другой шаг кроме 1 получим ошибку если далее будет TEST вместо CMP, с CMP же шаг может быть любым
        cmp     ESI,ECX
        jc      @sum2

Т.е. часто удобство программиста важнее скорости.

-- 02.12.2017, 16:26 --

(Оффтоп)

PS. Только у меня одного при помещении ASM кода в тег off пропадает кнопочка Развернуть - «Верните кнопку Развернуть»? А, похоже не только, уже упоминали - «Конфликт тегов off и syntax?».