Перемножение матриц - конкурс для программистов

[venco]

Zealint в сообщении #291525 писал(а):

Хорошо, может я чего-то не понимаю...

Если честно, я тоже не понимаю, что хочет сказать Dongara. То он говорит, что оптимальный алгоритм может умножить за 3 секунды, а всё, что было получено на конкурсе - фигня. То он возмущается тем, что такого не может быть и в 7.5 секунд уложиться невозможно. Это как я смог понять его туманные намёки. Будет лучше, если он прямо напишет, что он имеет в виду, и желательно с учётом поставленной задачи - умножение 16-ти битных целых чисел, а не вещественных.

[venco]

Вот код, который занимает большую часть времени.
Для матрицы 5000х5000 он вызывается $\frac {7^4\cdot 320^2} 4$ раз с $N=320$.
В среднем один цикл занимает 6 тактов и выполняет 32 умножения+сложения.
Буду рад, если кто-нибудь сможет его заметно ускорить.

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

#include <emmintrin.h>

// Calls instruction PHADDD from SSSE3
// For unknown reason it's not defined in the emmintrin.h
inline __m128i _mm_hadd_epi32 (__m128i a, __m128i b)
{
    __asm("phaddd %[b], %[a]" : [a]"=&x"(a): "[a]"(a), [b]"x"(b));
    return a;
}

// Multiply row of N short values pointed by a
// by 4 rows of N short values pointed by b.
// Difference between the rows is specified by bsize (in bytes).
// So, the b rows' pointers are:
//   b
//   (const short*)((const char*)b+bsize)
//   (const short*)((const char*)b+2*bsize)
//   (const short*)((const char*)b+3*bsize)
//
// Products for all 4 rows are then accumulated as int values and stored
// into 4 sequential int values pointed by m.
// N>0, and is multiple of 8.
// The pointers m, a, b, and the offset bsize are 16-byte aligned.

inline void row4_mul(int N, int* m, const short* a, const short* b, int bsize)
{
    // s0, s1, s2, s3 are accumulators for the products
    __m128i s0 = _mm_setzero_si128(), s1 = _mm_setzero_si128();
    __m128i s2 = _mm_setzero_si128(), s3 = _mm_setzero_si128();
    // v0, v1, v2, v3 are temporary registers
    __m128i v0, v1, v2, v3;
    // prepare pointers and offsets to 
    const short* be = b + N; // end pointer
    int b2a = (const char*)a - (const char*)b; // offset of a relative to b
    int bsize3 = 3*bsize; // offset of 3rd row in b
    do {
        // multiple and add 8 input values from each of the 4 rows at a time
        asm("movdqa  (%[b],%[b2a],1),%[v3] \n\t"
            "movdqa  (%[b]),%[v0] \n\t"
            "pmaddwd %[v3],%[v0] \n\t"
            "movdqa  (%[b],%[bsize],1),%[v1] \n\t"
            "pmaddwd %[v3],%[v1] \n\t"
            "movdqa  (%[b],%[bsize],2),%[v2] \n\t"
            "pmaddwd %[v3],%[v2] \n\t"
            "pmaddwd (%[b],%[bsize3],1),%[v3] \n\t"
            "add     $0x10,%[b] \n\t"
            "paddd   %[v0],%[s0] \n\t"
            "paddd   %[v1],%[s1] \n\t"
            "paddd   %[v2],%[s2] \n\t"
            "paddd   %[v3],%[s3] \n\t"
            : [s0]"=&x"(s0), [s1]"=&x"(s1), [s2]"=&x"(s2), [s3]"=&x"(s3),
              [v0]"=&x"(v0), [v1]"=&x"(v1), [v2]"=&x"(v2), [v3]"=&x"(v3),
              [b]"=&r"(b)
            : "[s0]"(s0), "[s1]"(s1), "[s2]"(s2), "[s3]"(s3),
              [b2a]"r"(b2a), [bsize]"r"(bsize), [bsize3]"r"(bsize3));
    } while ( b < be );
    // add accumulators into a single xmm register
    s0 = _mm_hadd_epi32(s0, s1);
    s2 = _mm_hadd_epi32(s2, s3);
    s0 = _mm_hadd_epi32(s0, s2);
    _mm_store_si128((__m128i*)m, s0);
}

void foo(int N, int* m, const short* a, const short* b, int bsize)
{
    row4_mul(N, m, a, b, bsize);
}
 

[Zealint]

venco в сообщении #291606 писал(а):

Вот код, который занимает большую часть времени.
Для матрицы 5000х5000 он вызывается $\frac {7^4\cdot 320^2} 4$ раз с $N=320$.
В среднем один цикл занимает 6 тактов и выполняет 32 умножения+сложения.
Буду рад, если кто-нибудь сможет его заметно ускорить.

У Вас получается 1.25 обращений к памяти на каждые 8 умножений, так что узким местом может оказаться именно это. Почему Вы уверены, что тут 6 тактов? Мне кажется, что ускорять надо не этот кусок, а его вместе с предыдущем циклом. Переорганизовав порядок умножения как у меня, получается 0.75 обращений на 8 умножений, но у меня такой подход оказался медленнее. И более того, такая же процедура (с точностью до порядка следования команд) была и у меня, и работала просто ужасно. Причина большой скорости у Вас явно в другом месте. Как-то Вы аккуратнее с кэшем работаете что-ли. Я, конечно же, разберусь.

[venco]

Zealint в сообщении #291710 писал(а):

Почему Вы уверены, что тут 6 тактов?

Определил из времени выполнения.

Zealint в сообщении #291710 писал(а):

Переорганизовав порядок умножения как у меня, получается 0.75 обращений на 8 умножений, но у меня такой подход оказался медленнее.

Интересно, каким образом. У меня либо не хватает регистров, либо приходится внутри цикла вместо paddd использовать более медленное суммирование через phaddd.

Zealint в сообщении #291710 писал(а):

Причина большой скорости у Вас явно в другом месте. Как-то Вы аккуратнее с кэшем работаете что-ли.

У меня эта функция вызывается так, что массивы b почти всегда находятся в кеше первого уровня. Из памяти (а точнее из кеша 2-го уровня) читается только массив a.

[Zealint]

venco в сообщении #291713 писал(а):

Определил из времени выполнения.

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

Цитата:

Интересно, каким образом. У меня либо не хватает регистров, либо приходится внутри цикла вместо paddd использовать более медленное суммирование через phaddd.

Именно. С регистрами надо колдовать, чтобы их хватило. Я написал в блоге, что именно происходит в моей процедуре: большие блоки делятся на блоки 2x2, которые перемножаются за 1 такт (8 умножений + 4 сложения через paddwd). За одно обращение к памяти извлекаются сразу 2 матрицы 2x2, за второе - еще две, и за третье - еще две. Далее, первые две умножаются на вторые две, а потом на третьи две, итого 32 умножения и всего три обращения. Всего 0.75 обращений на 8 умножений. Пришлось подумать над экономией регистров, по-моему даже пришлось сделать один лишний movdqa между регистрами, ради этой экономии. У меня подозрение, что моя кубическая версия будет быстрее, если я все-таки отображение файла в память сделаю. Надо проверить...

Цитата:

У меня эта функция вызывается так, что массивы b почти всегда находятся в кеше первого уровня. Из памяти (а точнее из кеша 2-го уровня) читается только массив a.

Кэш L2 находится вплотную к процессору, поэтому это еще не память в прямом смысле. Но в целом тогда мне понятно, почему 1.25 обращений к памяти у Вас быстрее. Интересно получается, надо объединить мой подход и Ваш, может ещё можно кубическую версию из 10 секунд вытащить. Я подумаю над этим чуть позже.

[venco]

Zealint в сообщении #291714 писал(а):

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

Посмотрел ваш код. Вы ничего не предприняли для улучшения работы кеша. Вам помогло то, что матрицы 320х320, которые реально умножаются после Штрассена, помещаются в кеш 2-го уровня. В моём коде порядок умножения такой, что кеш 1-го уровня активно используется.
А именно, я умножаю сначала первые 32 строки матрицы B на все строки матрицы A, потом следующие 32 строки матрицы B, и т.д. В результате практически всё время нужные элементы матрицы B лежат в кеше 1-го уровня, а не 2-го, как у вас. Ещё лучше этот метод работает, если матрица не помещается целиком даже в кеш 2-го уровня, как если не использовать Штрассен.

Zealint в сообщении #291714 писал(а):

Именно. С регистрами надо колдовать, чтобы их хватило. Я написал в блоге, что именно происходит в моей процедуре: большие блоки делятся на блоки 2x2, которые перемножаются за 1 такт (8 умножений + 4 сложения через paddwd). За одно обращение к памяти извлекаются сразу 2 матрицы 2x2, за второе - еще две, и за третье - еще две. Далее, первые две умножаются на вторые две, а потом на третьи две, итого 32 умножения и всего три обращения. Всего 0.75 обращений на 8 умножений.

В результате ваш цикл разбавлен кучей инструкций pshufd, а они тоже занимают время.
Сравните, на каждые 4 полезные pmaddwd и 4 полезные paddd лишних операций:
в моём коде: 5 чтений из памяти (кеша), причём одно чтение совмещено с pmaddwd,
в вашем: 3 чтения, 6 pshufd, и 2 пересылки в xmm регистрах.

[Zealint]

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

Добавлю: хотя моя кубическая версия неэффективна при работе с кэшем, она всего на секунду уступала Вашей (у меня было 14,5), причем был медленный ввод / вывод данных. Это даёт надежду на то, что, учтя Ваши замечания, можно написать довольно эффективную программу.

[venco]

Zealint в сообщении #291893 писал(а):

Я согласен с замечаниями, более того, приму меры для исправления ситуации. Кстати, эксперименты показали, что все мои "лишние" инструкции работают гораздо быстрее, чем подобная реализация с умножением по строкам.

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

Zealint в сообщении #291893 писал(а):

Добавлю: хотя моя кубическая версия неэффективна при работе с кэшем, она всего на секунду уступала Вашей (у меня было 14,5), причем был медленный ввод / вывод данных.

В моём "кубическом" решении умножение было написано без ассемблера - на GCC intrinsic functions (_mm_padd_epi16() etc). В конце конкурса я заметил, что компилятор напихал в цикл 4 лишние пересылки между xmm регистрами. Когда я их убрал, чистое умножение ускорилось почти в полтора раза.

[Dongara]

Zealint писал(а):

Правильно ли я понял, что Ваша программа не ускоряется методом Штрассена?

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

P.S.
Зря я ввязался в эту бучу: не надо было этого делать. Уже давно заметил: когда покидаешь привычный круг общения, то это приводит только к одним неприятностям.

[venco]

Dongara в сообщении #292006 писал(а):

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

Вы тоже участвовали в конкурсе?

[Zealint]

Dongara в сообщении #292006 писал(а):

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

Тов. Dongara, правильно ли я понял, что у Вас есть эффективный алгоритм, но Вы не хотите его "палить", боясь, что кто-то украдет идею? То есть Вы, видимо, хотите свой алгоритм кому-то продать (например, тем, кто делал BLAS)? Тогда прошу объяснить: что неправильного было сделано в конкурсе, раз Вы так накинулись на меня (и некоторых других) с пустыми обвинениями?

Dongara в сообщении #292006 писал(а):

P.S.
Зря я ввязался в эту бучу: не надо было этого делать. Уже давно заметил: когда покидаешь привычный круг общения, то это приводит только к одним неприятностям.

Совершенно верно, если нечего сказать по-существу, то лучше промолчать, а иначе получается именно то, что Вы имеете возможность наблюдать. Вот например, я видел, как люди соревнуются по одной задаче, в которой я разбираюсь (пока не хочу говорить, в какой). Их программы работали в несколько ТЫСЯЧ раз медленее моей. Но я даже не подал виду, что знаком с этой областью очень хорошо. Просто я тоже не хочу "палить" идею, пока не доведу её до строгого обоснования. А когда доведу мне пофиг, кто её украдет и кому продаст. Мне это безразлично.

[Dongara]

Zealint в сообщении #292012 писал(а):

видимо, хотите свой алгоритм кому-то продать

Целочисленный алгоритм без ограничений на элементы матрицы писался под заказ.

Zealint в сообщении #292012 писал(а):

Dongara в сообщении #292006 писал(а):

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

Совершенно верно.

Это я в основном про Вас: вот один известный человек поспорил с не менее известным товарищем по поводу алгоритма перемножения матриц для 65 нм. процессора под x32 (там одна нужная команда еще к вещественному порту была привязана). Один из них выиграл, но у него на это ушел примерно год. И это профессионалы. И если Вам нужен соответствующий алгоритм профессионального уровня, то надо и искать соответствующих специалистов, а не объявлять клоунские конкурсы: настоящие спецы над ним смеются, но, как говорят, "в тряпочку". И мне говорили: не ввязывайся.

[Zealint]

Dongara в сообщении #292054 писал(а):

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

Я поздравляю Ваших профессионалов и не сомневаюсь в их мастерстве.
Теперь по существу: я наконец-то Вас понял, но Вам придется поверить (хотя мне дела до этого нет), что конкурс этот придуман не зря, и мотивов Вам моих не понять, хотя намеков я сделал уже столько... в общем много. Вам правильно сказали, что ввязываться не нужно. Профессионалы, такие как проф. Грановский, скажем, мне не нужны. Далее, меня вообще мало интересуют простые задачи. Еще далее, меня не интересует оптимизация в 2, 5, 10, 100 раз и так далее (это ерунда), я занимаюсь задачами, которые нужно ускорять в миллионы раз; а зачем мне нужны были матрицы - это совершенно не важно. Знаете, специалисты в той области, которой я сейчас занимаюсь тоже без улыбки на Ваши ускорения алгоритмов из MKL всего на 50% (да хоть на 50 000%) смотреть не могут (не буду цитировать их слова). Именно поэтому я Вам сразу (найдите в этой теме) сказал: не вмешивайтесь в то, чего не понимаете (я-то сообразил, что больную для Вас тему поднял). Меня тоже предупредили, что кое-кто будет наезжать. Но я результатами конкурса доволен. И дело не совсем в цифрах. Думаю, разговор теперь окончен?

PS. Если все-таки ещё не поняли, то давайте смотреть друг на друга как на детей, которые просто развлекаются. Т. е. Ваши работы для меня - детский сад. Мои действия носят аналогичных характер для Вас. Добро?

[venco]

Dongara, а давайте вы будете за свои слова отвечать. Вот вопросы, которые вы проигнорировали:
http://dxdy.ru/post289575.html#p289575
http://dxdy.ru/post291532.html#p291532
http://dxdy.ru/post291567.html#p291567
http://dxdy.ru/post292010.html#p292010

[Dongara]

К тому, что я уже написал, мне добавить нечего. Успехов.