Блочное умножение матриц

Продолжение.
Этот пост - продолжение предыдущего, где обсуждалась производительность кода, генерируемого LLVM, по сравнению с GCC и ICL на примере умножения матриц. Тестировались реализации без блокирования. Здесь представлены результаты измерения производительности блокирующего алгоритма и пример автоматической генерации кода.
Поблочное умножение.
Алгоритм реализован для произвольных размеров матриц и блока. Производительность сравнивается с реализаций на asm как самой быстрой.

Ничего удивительного на графике не найти. Как бы ни был быстр asm код, обработка больших матриц и доступ к данным за пределами кэша ведет к замедлению. Видно падение производительности для объемов данных, превышающих 4 MB. Обработка по блокам, в основном, имеет дело с данными, лежащими в кэше, поэтому быстрее и, как видим, значительно. Умножение блока выполняется тем же самым asm кодом, с которым сравниваем. Скачки для блочной реализации объясняются разными размерами блока. При произвольности параметров есть условия максимальной производительности:

• размеры матриц кратны размеру блока,
• размер блока кратен 8-ми элементам (два SSE2),
• указатели выравнены на 16 байтов.

Использовал размер блока до 256х256. Каждый раз максимальная производительность достигалась при наибольшем размере блока. Небольшой размер блока дает гибкость в обработке (бо'льшая часть матриц может быть накрыта блоками), но добавляет больше непроизводительных операций (копирование данных в/из блока)
Код можно генерировать
Может оказаться интересным автоматическая генерация кода, возможная, например, при использовании metaocaml. Сгенерировать код умножения матриц с разверткой циклов обнуления массива результата и внутреннего цикла умножения элементов матриц можно с помощью такого кода

Использованы 2 функции развертки: одна разворачивает весь цикл, другая разворачивает цикл до 256 с фактором. Развертка, как видим, выполняется на 8. Размер блока фиксирован - 256х256. Мета-средства, а именно, многостадийное программирование multistage programming , нужны здесь для генерации тела операции, которая размножается. В коде таких две: присвоение нулю и умножение элементов. С помощью офшоринг-механизма ".!{Trx.run_gcc}" генерируется следующий С-код

void mmm_256_u(float*a_4,float*b_5,float*c_6)
{
 int i_7 ;
 for(i_7 = 0; i_7 <= 255; i_7++)
 {
   int i256_8 ;
   int ii_13 ;
   int k_9 ;
   i256_8 = i_7 * 256;
   for(ii_13 = 0; ii_13 <= 32 - 1; ii_13++)
   {
     c_6[i256_8 + (ii_13 * 8 + 0)] = 0.;
     c_6[i256_8 + (ii_13 * 8 + 1)] = 0.;
     c_6[i256_8 + (ii_13 * 8 + 2)] = 0.;
     c_6[i256_8 + (ii_13 * 8 + 3)] = 0.;
     c_6[i256_8 + (ii_13 * 8 + 4)] = 0.;
     c_6[i256_8 + (ii_13 * 8 + 5)] = 0.;
     c_6[i256_8 + (ii_13 * 8 + 6)] = 0.;
     c_6[i256_8 + (ii_13 * 8 + 7)] = 0.;
   } ;
   for(k_9 = 0; k_9 <= 255; k_9++)
   {
     int k256_10 ;
     float aa_11 ;
     int ii_12 ;
     k256_10 = k_9 * 256;
     aa_11 = a_4[i256_8 + k_9];
     for(ii_12 = 0; ii_12 <= 32 - 1; ii_12++)
     {
c_6[i256_8 + (ii_12 * 8 + 0)] = c_6[i256_8 + (ii_12 * 8 + 0)] + aa_11 * b_5[k256_10 + (ii_12 * 8 + 0)];
c_6[i256_8 + (ii_12 * 8 + 1)] = c_6[i256_8 + (ii_12 * 8 + 1)] + aa_11 * b_5[k256_10 + (ii_12 * 8 + 1)];
c_6[i256_8 + (ii_12 * 8 + 2)] = c_6[i256_8 + (ii_12 * 8 + 2)] + aa_11 * b_5[k256_10 + (ii_12 * 8 + 2)];
c_6[i256_8 + (ii_12 * 8 + 3)] = c_6[i256_8 + (ii_12 * 8 + 3)] + aa_11 * b_5[k256_10 + (ii_12 * 8 + 3)];
c_6[i256_8 + (ii_12 * 8 + 4)] = c_6[i256_8 + (ii_12 * 8 + 4)] + aa_11 * b_5[k256_10 + (ii_12 * 8 + 4)];
c_6[i256_8 + (ii_12 * 8 + 5)] = c_6[i256_8 + (ii_12 * 8 + 5)] + aa_11 * b_5[k256_10 + (ii_12 * 8 + 5)];
c_6[i256_8 + (ii_12 * 8 + 6)] = c_6[i256_8 + (ii_12 * 8 + 6)] + aa_11 * b_5[k256_10 + (ii_12 * 8 + 6)];
c_6[i256_8 + (ii_12 * 8 + 7)] = c_6[i256_8 + (ii_12 * 8 + 7)] + aa_11 * b_5[k256_10 + (ii_12 * 8 + 7)];
     } ;
   } ;
 } ;
}

Пришлось отредактировать только тип данных, заменил double на float - ocaml знает только тип данных с плавающей запятой double. Производительность кода с разверткой плохая. Как мы уже упоминали, ручная развертка мешает Intel компилятору векторизовать код. Следовательно, средства мета-программирования, использованные здесь для развертки, оказались неприменимы.

Высокой производительности достигает код без развертки, генерируемый из такого ocaml кода

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

для матриц размером 2072x2072 и блоков размером 256x256

• unrolled: 1161 MFlops,
• vectorized by ICL: 2339 MFlops.

Многостадийное программирование особенно эффективно в тех случаях, когда на этапе компиляции известна дополнительная информация о данных. Например, есть умножения на 1.0 или 0.0, или вычисляются значения тригонометрических функций для pi/2, pi/4. Примером такой задачи является FFT, решение ее подробно рассмотрено в работе "A Methodology for Generating Verified Combinatorial Circuits" авторами O.Kiselyov, K.N.Swadi и W.Taha - восхитительный образец использования мета-программирования.