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

[Dmitriy40]

Может и не всегда, однако в большинстве случаев. И проигрыш от 30% до 12% (в умножении) несравним с проигрышем в 10-15 раз - и это после неслабых усилий по оптимизации.
Потому, если надо "ехать", то выгоднее применить библиотеку, которая почти гарантированно выдаст далеко не худший результат в любых условиях, чем писать что-то своё и долго потом оптимизировать и адаптировать под разные внешние условия. Пусть даже результат библиотеки будет и не 100% оптимальным - ради этого оптимума придётся потратить месяцы (а то и годы) труда. Имхо.

[YYSS]

Dmitriy40
Да я не против: просто хотел подчеркнуть, что планка высших достижений иногда неподвластна даже очень сильным коллективам, которым безусловно является Intel, теоретическую поддержку которому оказывает нехилый профессорский коллектив ученых. И потом, для алгоритмов диагонализации упакованных матриц разность достигает более 50-ти процентов. И не надо хватать звезды с неба, чтобы реализовать довольно примитивный алгоритм перемножения матриц, уступающий в скорости dgemm Intel MKL всего процентов 30.
Что касается Relatively Robust Representations алгоритма, то разрыв достигает 2000%!

[Dmitriy40]

YYSS
Да я тоже не против.
К тому же здесь речь шла лишь об умножении матриц.
Кстати у меня вот пока (со вчерашнего дня) не получилось догнать MKL даже близко, никак не придумаю как "кормить" fmadd команды с темпом 16 double чисел за такт и при этом не проседать из-за ассоциативности кэша. Асимптотика показывает (уже раз 100 пересчитывал в разных вариантах) что это реально, надо лишь правильно сделать разбиение задачи. Так что правы когда говорят что это хороший учебный пример на многостадийную оптимизацию кода, я уже насчитал полдесятка методов ускорения, а до пиковой скорости ещё далеко. Так что насчёт звёзд с неба Вы чуток погорячились думаю, не так всё просто. Подсматривать же в чужом коде неспортивно/неинтересно.
На этой оптимистичной ноте можно пожалуй и закончить отступление.

[Geen]

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

библиотеку Intel MKL даже близко не догнали по скорости

Если я правильно помню, эта библиотека оптимизирована под архитектурные особенности одноимённых процессоров...

[Dmitriy40]

Geen
Я не уверен точно т.к. не смотрел её, но из очевидных соображений думаю она или сама подстраивается в момент выполнения под текущий процессор, или тип процессора задаётся ключами компилятору (вероятней всего). И да, разумеется оптимизирована. Ещё встречал утверждения что на процессорах AMD она вовсе не такая быстрая даже по сравнению с аналогами.

[YYSS]

Geen в сообщении #1225922 писал(а):

Если я правильно помню, эта библиотека оптимизирована под архитектурные особенности одноимённых процессоров...

Совершенно верно: она определяет тип процессора (интеловского!) и подставляет соответствующий код: именно поэтому dgemm Intel MKL таких агромадных размеров, исчисляющихся мегабайтами.

[bondkim137]

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

не проседать из-за ассоциативности кэша

Может, все-таки всю матрицу стоит повернуть ? Или хотя бы пачку столбцов так, что б шириной с линейку кеша получилось.

[Dmitriy40]

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

Может, все-таки всю матрицу стоит повернуть ? Или хотя бы пачку столбцов так, что б шириной с линейку кеша получилось.

Это разные практически независимые оптимизации. И поворот всей матрицы или только текущего столбца (группы столбцов) эквиваленты.
Чтобы снизить трафик памяти надо в кэше перемножать кусочки матриц не менее 128х128, а ассоциативность кэша всего 8, поэтому вероятно надо перемножать неквадратный блок, типа 8х2048, что может быть больше самой матрицы ... В общем хорошая задачка на организацию потоков данных с учётом возможностей аппаратуры.

[Pavia]

Я бы тоже транспонировал.
Есть быстрый код, для 4х4
http://masm32.com/board/index.php?PHPSE ... 6#msg65026
И как пишут для любого размера.
http://masm32.com/board/index.php?topic=6140.0

[YYSS]

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

YYSS
К тому же здесь речь шла лишь об умножении матриц.

Дело в том, что умножении матриц играет важную вспомогательную роль во многих ключевых алгоритмах линейной алгебры, например, в диагонализации, интернетные примеры для которой я и привел. Т.е. нет скорости умножения матриц - нет скорости и соответствующего алгоритма. Хотя в моих ссылках можно еще много чего интересного найти, касающегося проколов Intel MKL в их, как они рекламируют, самых быстрых алгоритмах.

Pavia
Есть такая программа: PC GAMESS. В ней есть ссылка на ассемблерный код, который содержит и быстрое перемножение матриц, не уступающее по скорости dgemm Intel MKL.

[Pavia]

YYSS
Я тут посмотрел http://eigen.tuxfamily.org/ она на Си++ но оптимизирована и под SSE, AVX, CUDA и NEON.
Там и перемножение и транспонирование есть, не надо голову ломать.

[Andrey_Kireew]

Для интереса привожу код, который работает быстрее чем с транспонированием,
правда не намного, но тем не менее (800 мс против 1100 мс)

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

int jj,kk,i;
   int nb=100;
   
   for (jj=0; jj<NStr; jj=jj+nb)
   for (ii=0; ii<NStr; ii=ii+nb)
   for (j=jj; j<jj+nb; j++)
   for (i=ii; i<ii+nb; i++)
       *(D+i+j*NStr)=0;
   
   for (kk=0; kk<NStr; kk=kk+nb)
   for (jj=0; jj<NStr; jj=jj+nb)
   for (ii=0; ii<NStr; ii=ii+nb)
       
   for (k=kk; k<kk+nb; k++)    
   for (j=jj; j<jj+nb; j++)
   for (i=ii; i<ii+nb; i++)
       
        *(D+i+j*NStr)+=*(B+i+k*NStr)**(M+k+j*NStr);
 

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

[bondkim137]

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

[Andrey_Kireew]

Там не менее следующий код работает в 1.5 раза медленнее предыдущего (и это не считая накладных расходов на транспонирование)

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

 int jj,kk,i;
   int nb=100;
   
   for (jj=0; jj<NStr; jj=jj+nb)
   for (ii=0; ii<NStr; ii=ii+nb)
   for (j=jj; j<jj+nb; j++)
   for (i=ii; i<ii+nb; i++)
       *(D+i+j*NStr)=0;
   
   
   
   for (ii=0; ii<NStr; ii=ii+nb)
   for (jj=0; jj<NStr; jj=jj+nb)    
   for (kk=0; kk<NStr; kk=kk+nb)
       
   for (i=ii; i<ii+nb; i++)    
   for (j=jj; j<jj+nb; j++)
   { val=0;
   for (k=kk; k<kk+nb; k++) 
       val+=*(BT+k+i*NStr)**(M+k+j*NStr);
       *(D+i+j*NStr)=val;
   }
 

так, что предыдущий код самый быстрый из того что получилось

[YYSS]

Andrey_Kireew

Советую вам ознакомиться с основами оптимизации циклов здесь.