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

Здравствуйте товарищи-алгоритмисты. Хочу поделится тут с вами своими (и не только) соображениями по поводу матричного умножения. Тема эта как казалось мне уже давно избитая, но когда начинаешь копать глубже при оптимизации умножения всплывают всё новые и новые факты. Вообщем стоит задача: написать высокопроизводительное матричное умножение так чтобы вычисления с плав. точкой производились на пиковой (или почти пиковой) скорости процессора (в расчет не берется многоядерность). Матрицы берутся общего вида, достаточно большого размера (не помещаются в L2 кэше данных). Т.е. в матричном виде имеем:
C = C + A*B, где A,B,C есть M x K, K x N и M x N матрицы.
Высокая производительность такого умножения вообще достигается двумя путями:
1) Оптимизация внутреннего цикла, за счет использования векторизации и блокирования регистров – т.е. сокращения кол-ва перемещений данных между кэшем и регистрами. Последнее достигается за счет использования всех доступных xmm – регистров (8 под x86 и 16 под x64) и разворотом внутреннего цикла. При этом данные, загруженные в регистры, повторно используются. Обычно половина из доступных регистров используется для накапливания сумм произведений, остальные для вычислений. Этот вопрос уже поднимался мной в этой теме и впринципе ответ на него я получил.
2) Блокирование кэша – сокращение IO – трафика между ОЗУ и кэшем. Идея заключается в том, чтобы загружать блоки матриц в кэш и там эффективно их перемножать, при этом чтобы данные из одного или двух блоков максимально долго оставались в кэше и многократно использовались. Главный вопрос, из-за которого я и создал эту тему, состоит в том: как выбирать оптимальные размеры блоков и саму схему блокирования кэша, так чтобы вычисления производились на пике вычислений с плав. точкой? При том, такая схема должна эффективно работать на большинстве существующих архитектур процессоров и не должна зависеть от какого–то конкретного типа процессора.
Вообще есть ещё и 3-ий путь – это уменьшение кол-ва умножений за счет увеличения кол-ва сложений (алгоритмы типа Штрассена). Но здесь мы рассматривать их не будем так как это отдельная тема для исследований и этот путь идет в разрез с путем блокирования кэша. Хотелось бы достичь высокой производительности алгоритма имеено за счет первых двух. И это можно сделать. В подтверждение тому выскопроизводительные реализации GEMM из оптимизированных BLAS-библиотек (INTEL MKL, ATLAS и GotoBLAS).
Теперь рассмотрим более подробно 2-ой вопрос.
Первое что приходит на ум, это взять 3 блока A,B,C примерно одной размерности, загрузить их в кэш и там перемножить. При этом размер блока берется таким, чтобы эти три блока помещались в L2 кэш. (Lb ≤ sqrt(L2/3))

Тогда алгоритм блочного умножения матриц будет выглядеть примерно так:

Nb = N/Lb;
For (ib = 0;ib<Nb;ib++)
  For (jb = 0;jb<Nb;jb++)
    For (kb = 0;kb<Nb;kb++)
         Cibjb += Aibkb*Bkbjb; перемножение блоков

При стандартной схеме (“строка - столбец”) каждое перемножение первой строки из блока матрицы A на столбцы из блока матрицы B будет самым дорогостоящим, так как в это время ни элементов из блока A, ни элементов из блока B ещё нет в кэше. Умножение 2-ой и последующих строк блока A на столбцы из того же блока B будет происходить уже гораздо быстрее, так как загружаться в кэш из основной памяти будут только строки блока A а весь блок B уже будет находится полностью в кэше.
При такой схеме наиболее долго в кэше будет оставаться блок матрицы C, а элементы из блока матрицы A будут использоваться лишь единожды. Т.е. здесь блокирование матрицы A не имеет смысла.

Также очень важным для повышения производительности алгоритма является то, чтобы доступ к данным в памяти был строго последовательным (чтобы элементы матриц читались из строк а не из столбцов).
В случае когда матрицы квадратные можно применить транспонирование матрицы B и перемножать схемой “строка на строку”. Такая схема годится для квадратных матриц одного размера. Но когда матрицы общего вида, применять транспонирование матрицы B уже нельзя. Чтобы доступ к элементам был всё время последовательным можно применить алгоритм перемножения блоков матриц, переставив j-ый и k-ый циклы местами:

For (i = 0;i<Lb;ib++)
  For (k = 0;k<Lb;k++)
    For ( j = 0;j<Lb;j++)
        cij += aik*bkj;

Этот вариант уже лучше (в плане производительности) но не получается достичь пиковой скорости на одном ядре. Даже с оптимизацией внутреннего цикла (векторизация и блокирование регистров) достигается только 50% от пика.
Скорее всего квадратно-блочная схема блокирования далеко не оптимальная. В случае квадратных матриц, был также рассмотрен другой вариант блокирования, предложенный товарищем Leo. Это вариант, когда блоки берутся не квадратные, а строчные во всю длину строк матриц. При этом большой строчный блок транспонированной матрицы B заполняет значительную часть L2 кэша и удерживается там достаточно долго, а строки матрицы A меняются одна за другой . После того как вся матрица A будет помножена на блок B, в L2 кэш загружается следующий строчный блок матрицы B и всё повторяется.

Алгоритм выглядит так:

For (ib=0;ib<Nb;ib++)
  For (j=0;j<N;j++)
       Cj,ib += Aj * Bib // умножение строки матрицы A на блок матрицы B

Но и этот вариант также далек от оптимального. Этот алгоритм также выжимает только половину от пиковой производительности процессора.
Далее я решил узнать как осуществляется блокирование кэша в оптимизированных BLAS – библиотеках и наткнулся на одну очень интересную статью товарища Goto, автора высопроизводительной библиотеки GotoBlas. Статья называется “Anatomy of High-Performance Matrix Multiplication”. В ней есть все ответы на мои вопросы: там подробно описывается как осуществлять блокирование кэша и как выбирать размеры блоков так, чтобы достичь максимальной производительности при умножении матриц. Основная идея там заключается в том, чтобы разложить GEMM на вложенные “уровни ”. Каждый “уровень” это какая –то процедура. Например, один из возможных вариантов выглядит так:

GEMM раскладывается на несколько вызовов GEPP(умножение столбцового блока A –“панель” на строчный блок – также “панель”). В свою очередь GEPP может быть разложена на несколько вызовов GEBP (умножение прямоугольного блока на строчный).
И если, достичь максимальной производительности от процедуры последнего уровня (от GEBP), то GEMM также будет выполнятся на максимальной скорости процессора.
Очень советую прочесть эту статью чтобы быть в курсе.(есть её собственный перевод на русский язык). И далее в статье делается упор на рассмотрение именно этой процедуры последнего уровня. Если рассмотреть GEBP более подробно то видно:

что прямоугольный блок (mc * kc) матрицы A (A1), и два nr – столбца из B (bi) и С(ci) должны находится в L2 кэше. Причем (mc * kc) матрица должна удерживаться там достаточно долгое время(пока будет нужна). Кроме того данные, помещенные в L2 кэш должны быть адресованы в TLB L2 кэш. И адреса (mc * kc) матрицы также должны удерживаться в L2 TLB кэше. Если 2-ое условие не будет соблюдаться , то будет происходить большее кол-во TLB – промахов. А TLB-промах, как пишет Goto, гораздо дороже кэш промаха. Т.е. имеем:

(1) mc * kc  + 2*(mc*nr + kc* nr) <=L2  
и
(2) (mc * kc  + 2*(mc*nr + kc* nr)) записей <=TLB L2

Причина появления двойки в формулах (1) и (2) Goto объясняет следующем: когда потребуется загрузить следующие два сi+1 и bi+1 блочных столбца в L2 кэш и адресовать их L2 TLB, то должны быть заменены именно ci и bi столбцы и соответствующие записи в TLB. Однако, после завершения вычисления ci = ci + A1*bi может оказаться так что самые старые линейки и записи A1 будут заменены на новые bi+1 и ci+1. Коэффициент “2” позволяет хранить линейки и записи, связанные с bi и ci, вместе с линейками и записями A1, bi+1 , ci+1 до того момента , пока не потребуется загрузка и адресация следующих bi+2 и ci+2 , которые заменят уже не нужные данные и адреса из bi и ci.
Обычно, для вычисления параметров мс и kc используется меньшее из (1) и (2).

Кроме того, как я уже писал, доступ к данным должен осуществляться строго последовательный. Для этого Goto предлагает делать упаковку блоков –копировать специальным образом( в соответствии со схемой блокирования) непоследовательные данные из блока в последовательный одномерный массив.
Тогда алгоритм GEMM будет выглядеть примерно так:

For (k=0;k<K;k+=kc) {
  // выполнить k-раз GEPP
       B1 = Pack Bk ; // запаковать блок kc x N матрицы B в посл. массив B1
    For (j=0;j<M;j+=mc) {
          A1 = Pack Ajk ; // запаковать блок mc x kc матрицы A в посл массив A1
          Cj += A1*B1; // выполнить GEBP  
    }
}

Рассмотрим теперь как умножаются запакованные блоки и как они вообще упаковываются:

Блок матрицы A упаковывается так, чтобы каждая mr * kc подматрица блока хранилась последовательно в памяти. Блок матрицы B упаковывает последовательно kc * nr подматрицы. Параметры mr и nr ограничены кол-вом доступных регистров. Обычно половина доступных xmm – регистров используется для хранения mr * nr подматрицы Cj.
Оставшиеся регистры используются для умножения mr – элементов A1 на nr –элементов B1.
Тогда алгоритм GEBP:

For (i=0;i<N;i+=nr) {
  For (j=0;j<mc;j+=mr)
      C1 = A1j * B1i  // умножение mr * kc блочной строки на kc* nr блочный столбец
 C += unpack C1 ; // распаковать С1 в С
}

Стоит сказать пару слов об произведении блочной строки на блочный столбец:
Из A1 в регистры последовательно загружается mr – элементов, Из B1 в регистры последовательно загружается nr – элементов. Далее они перемножаются и результат перемножения сохранятся в регистрах для накопления произведений. Далее на место старых в регистры загружаются следующие mr –элементы из A1. Точно также из B1 следующие nr – элементы. И всё повторяется. И так kc – раз.
После того как mr * kc блочная строка полностью умножена на kc* nr блочный столбец, сумма произведений, накопившаяся в регистрах, выгружается в последовательный массив C1. Таким образом, мы шагаем в памяти строго последовательно.
Из посл. картинки видно, что ещё одним условием для эффективной реализации GEBP умножения будет:

(3)  kc * nr + 2*(mr * nr + mr  * kc)  <= L1

Причина появления двойки такая же как и для условий (1) и (2) : удержание всего kc * nr блочного столбца в L1.
Зная кол-во доступных регистров можно определить mr и nr. (nr также определяется, исходя из пропускной способности между L2 и регистрами и скоростью вычислений с плав. точкой – см. более подробно в статье).
Для 8-ми доступных xmm регистров mr x nr обычно выбирается: 2x4, 2x2 и 4x2. Из условия(3) можно найти kc. А из условий (1-2) - mc.
Также Goto в своей статье сообщает, что на практике параметры, выбранные по условиям (1-3) не всегда оптимальны (на выбор этих параметров влияет ещё набор ассоциативности и стратегии замещения кэша) и предлагает выбирать их так:
kc –чисел с плав. точкой должны занять половину размера страницы памяти.
mc выбирается так, чтобы mc*kc матрица A1 занимала примерно половину меньшего из :
- памяти, адресуемой TLB L2;
- L2 кэша.
Например, для моего процессора AMD Athlon64 x2 (L1Cache = 64Kb, L2Cache = 512Kb, TLB L2 = 512 записей, pagesize = 4KB) я вычисляю их так:

mr и nr я взял 2x4.
kc = pagesize / 2 = 4096 Bytes / 2*8 Bytes = 256.
 
Проверка:
kc*nr + 2*(mr*nr +mr*kc)  <= L1
kc <= (L1-2*mr * nr) / (2*mr + nr) = (64 * 1024 – 2*2 * 4 * 8) Bytes / (2*2 +4) * 8 Bytes <=1022
256 <= 1022.
 
(1) (mc x kc) <= L2/2
mc <= L2 / 2 * kc = 512 * 1024 Bytes / 2 * 256 * 8 Bytes = 128
 
 (2) (mc x kc) <=TLB L2 /2
mc <= TLB L2/2*kc = 512*4096 Bytes / 2*256*8 Bytes = 512

Меньшее из (1) и (2) есть 128.

Таким образом для моего процессора у меня получились такие вот параметры:
mr = 2, nr = 4, kc = 256, mc = 128

Я запрограммировал этот алгоритм умножения матриц, предлагаемый Goto, но добиться максимальной производительности от него у меня опять не получилось. Снова те же 50% от пика. Также я пробовал брать mr x nr = 2x2, производительность аналогичная.
Хотя этот подход, предложенный Goto, позволяет достичь 90 и более % от максимальной производительности на большинстве архитектур. Это показывают тесты в статье да и сама библиотека GotoBlas даёт неплохие результаты, производительность которой сравнима с Intel MKL, а в некоторых случаях и превосходит её. Возможно, я что то упустил при реализации. Кому интересна эта тема советую прочитать статью и попробовать реализовать подход предложенный там, возможно у кого-то получится. Саму статью, её перевод и мою реализацию описанного в ней подхода можно скачать здесь.
Также было бы неплохо, если кто-нибудь предложит свой вариант блокирования и выбора размеров блоков.

Сообщение отредактировано: vitaly333 - 13.11.09, 21:56

Я бы использовал рекурсию,многопоточную реализацию и возможно модульную математику.Да и вот циклы For медленные,Я бы через DO сделал.
А вообше что по твоему "пик производительности" ты в диспечере задач штоли смотриш?.

Лично мне в последнее время нравится cache-oblivious подход к умножению матриц, в сочетании с использованием оптимизированных ядер для нижнего уровня рекурсии. Не надо танцев с бубном вокруг кэшей разных уровней и разного размера, учета таких тонкостей, как inclusive/non-inclusive кэш и т.д. Алгоритм сам использует кэш довольно-таки эффективным образом.

Цитата

А вообше что по твоему "пик производительности" ты в диспечере задач штоли смотриш

Она высчитывается в флопсах(операциях с плав. точкой в секунду):
ПП = [Частота 1 процессора]*[кол-во ядер]*[кол-во АЛУ]*[кол-во арифм. операций с плав. точкой, выполняемых за один такт процессора]

Цитата

Лично мне в последнее время нравится cache-oblivious подход к умножению матриц, в сочетании с использованием оптимизированных ядер для нижнего уровня рекурсии. Не надо танцев с бубном вокруг кэшей разных уровней и разного размера, учета таких тонкостей, как inclusive/non-inclusive кэш и т.д. Алгоритм сам использует кэш довольно-таки эффективным образом.

А поподробнее можно? Что за алгоритм? Какую производительность показывает?

Цитата vitaly333 @ 15.11.09, 21:37

Она высчитывается в флопсах(операциях с плав. точкой в секунду):
ПП = [Частота 1 процессора]*[кол-во ядер]*[кол-во АЛУ]*[кол-во арифм. операций с плав. точкой, выполняемых за один такт процессора]

Интересно, и какие же параметры (кроме частоты) ты берешь для своего атлона х2 ?
Надеюсь, понимаешь, что кол-во ядер имеет смысл учитывать только при многопоточном расчете, т.к. при однопоточном на x2 по этой формуле всегда будешь получать не более 50% ПП. С числом АЛУ и фп-операций за такт тоже не все так просто - нужно учитывать особенности архитектуры, в частности, как я уже упоминал SSE2 только в Core2 и AMD K10 выполняются за одну микрооперацию, а в более ранних моделях за две - если этого не учитывать, то также можно ошибиться ровно в 2 раза

В электронном виде можно найти 'Cache-Oblivious Algorithms', Harald Prokop. Там описан простейший рекурсивный алгоритм умножения матриц и другие алгоритмы. Правда, в оригинальном виде быстродействие у него очень низкое, но если не доводить рекурсию до конца (т.е. до умножения матриц 1х1), а спустившись до матриц 8х8 использовать эффективное ядро, то можно достичь неплохих результатов. 100% производительности процессора не выберешь, но зато полученный код будет работать с одной и той же эффективностью на любом процессоре - с любым размером кэша и без предварительной настройки.

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

Цитата andriano @ 16.11.09, 13:25

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

Нет, тут ещё два момента:
1. эффективная работа TLB (промахи TLB тоже оказывают существенное воздействие)
2. эффективная работа с уже загруженными в кэш матрицами, что тоже не совсем просто

Цитата

В электронном виде можно найти 'Cache-Oblivious Algorithms', Harald Prokop. Там описан простейший рекурсивный алгоритм умножения матриц и другие алгоритмы. Правда, в оригинальном виде быстродействие у него очень низкое, но если не доводить рекурсию до конца (т.е. до умножения матриц 1х1), а спустившись до матриц 8х8 использовать эффективное ядро, то можно достичь неплохих результатов. 100% производительности процессора не выберешь, но зато полученный код будет работать с одной и той же эффективностью на любом процессоре - с любым размером кэша и без предварительной настройки.

Посмотрю. Если вы делали его реализацию не могли бы вы привести какие-то конкретные цифры производительности, которую он показывает?

Цитата

Интересно, и какие же параметры (кроме частоты) ты берешь для своего атлона х2 ?

Для моего атлона я считал так(конечно многоядерность я не учитываю, т.к. речь идет о достижении пика на одном ядре):
ПП = [2412]*[1]*[1]*[2] = 4.8 Гфлопа
Реализации блочных алгоритмов показывают примерно 2 - 2.6 Гфлопа.
DGEMM из gotoBlas и MKL дают 3.8 - 4.4 Гфлопа (т.е. в любом случае, даже если я и ошибся с подсчетом, goto-алгоритм должен выдавать 80-90% а получается только около 50%).

Сообщение отредактировано: vitaly333 - 16.11.09, 18:32

Цитата vitaly333 @ 16.11.09, 18:27

Посмотрю. Если вы делали его реализацию не могли бы вы привести какие-то конкретные цифры производительности, которую он показывает?

Нет, я пока только планирую приступать к работе над ним...

Цитата vitaly333 @ 16.11.09, 18:27

Реализации блочных алгоритмов показывают примерно 2-2.6 Гфлопа.
DGEMM из gotoBlas и MKL дают 3.8 - 4.4 Гфлопа (т.е. в любом случае, даже если я и ошибся с подсчетом, goto-алгоритм должен выдавать 80-90% а получается только около 50%).

В таком случае ты мог ошибиться с реализацией goto-алгоритма

Цитата

В таком случае ты мог ошибиться с реализацией goto-алгоритма

Я вот про это и говорю, может где-то чего упустил... Много раз уже проверял, моя реализация сделана точно по статье... может быть товарищ Goto не упомянул ещё каких-нибудь деталей, влияющих на реализацию...незнаю сложно сказать

Покажи полностью прогу

Цитата vitaly333 @ 18.11.09, 22:49

Я вот про это и говорю, может где-то чего упустил... Много раз уже проверял, моя реализация сделана точно по статье...

В статье, например, не затрагиваются такие детали, как формат хранения матрицы. Судя по предыдущей теме, ты используешь double** с построчной аллокацией, а тот же albom - один непрерывный блок double*. Для самих больших матриц это не должно играть существенной роли (если только выравнивание строк не совпадет с алиасингом кэша), а вот для маленьких packed-матриц Goto это весьма существенно - они обязательно должны быть непрерывными едиными массивами, а не независимыми наборами строк как в double**

Сообщение отредактировано: leo - 19.11.09, 08:27

leo верно говориш,можно даже списком вроде они пошустрее матриц

Цитата

для маленьких packed-матриц Goto это весьма существенно - они обязательно должны быть непрерывными едиными массивами, а не независимыми наборами строк как в double**

Да так и есть, каждый блок запаковывается (копируется в соответствии с выбранной схемой блокирования регистров) из 2-мерного массива double** (исходной матрицы) в непрерывный одномерный рабочий массив double*.

Цитата

Покажи полностью прогу

Всё есть в ссылке которую я приводил выше.

Сообщение отредактировано: vitaly333 - 19.11.09, 12:52

Я реализовал программу(Visual C++ 2008) которая каждую матрицу разбивает на четыре блока и перемножает их между собой,как ее оптимизировать чтобы каждая матрицы разбивалась на оптимальное количество блоков,и производилось умножение?

Сообщение отредактировано: real27 - 13.03.10, 19:12

Прикреплённый файлishodnik.rar (0.55 Кбайт, скачиваний: 343)