Диссертация (Оптимизация размещения массивов в общей памяти), страница 9

Поэтому, если размерматрицы A не превосходит размера кеш-памяти ( N 2 * DataTypeSize  CacheSize ), токоличество кеш-промахов легко подсчитать и оно не превосходит числоN 2 * DataTypeSize. Поэтому будем предполагать, что размер матрицы A большеBlockSizeчем размер кеш-памяти.Подсчитаем количество кеш-промахов, возникающих при обращении квхождениям (1), (2), (3), (4) независимо:CacheMissCount = CacheMissCount(1) + CacheMissCount(2)+ CacheMissCount(3) + CacheMissCount(4)55На одной и той же итерации к вхождениям (2), (3), (4) при исполнениипрограммы происходят обращения раньше, чем вхождению (1). Из того, чтоиндексные выражения при вхождениях (1) и (2) совпадают, следует, чтоколичество промахов к (1) равно нулю:CacheMissCount(1) = 0Во вхождении (2) обращение к одним и тем же данным происходит, приизменении счетчика k.

Между двумя последовательными обращениями к одномуи тому же элементу матрицы A во вхождении (1) производит обращение ко всемэлементам матрицы A. Следовательно, при повторном обращении к элементу Aijбудет всегда происходить кеш-промах.CacheMissCount(2) N 3*DataTypeSizeBlockSizeКоличество кеш-промахов, порождаемых вхождением (3), не существенно,т.к. оно ограничено O(N 2 ) .Между двумя последовательными обращениями к одному и тому жеэлементу матрицы A во вхождении (4) происходит счиssтывание строки матрицыA. Подсчитаем, эффективное количество банков памяти (EffectiveBankCount) иэффективную ассоциативность кеш-памяти (EffectiveCacheAssoc) для этоговхождения.

В случае, когдаEffectiveCacheAssoc > CacheAssoc,выполняется равенствоCacheMissCount(4) =N2BlockSize, в противном случае выполняется равенствоCacheMissCount(4) =Таким образом:N.BlockSize56CacheMissCount N 3*DataTypeSizeBlockSize O(N 2 )Пример 8. Рассмотрим алгоритм Флойда-Уоршалла нахождения матрицыкратчайших расстояний, в котором матрица размещена блочно:Листинг 12. Алгоритм Флойда-Уоршалла с блочным размещением матрицыкратчайших расстояний.for (dk=0; dk<N; dk+=d)for (k=MAX(dk*d, 0); i < MIN((dk+1)*d, N); ++k)for (di=0; di<N; di+=d)for (di=0; di<N; di+=d)for (i=MAX(di*d, 0); i < MIN((di+1)*d, N); ++i)for (j=MAX(dj*d, 0); j < MIN((dj+1)*d, N); ++j)A[F2(di, dj, i, j)] (1) += MIN(A[F2(di, dj, i, j)] (2),A[F2(di, dk, i, k)] (3) + A[F2(dk, dj, k, j)] (4))Функция F2 принимает следующий вид:F2(di, dj,i,j)  di*d*N  d j*d 2  i*d  jКак и в предыдущем примере подсчитаем количество кеш-промахов,возникающих при обращении к вхождениям (1), (2), (3), (4) независимо:CacheMissCount = CacheMissCount(1) + CacheMissCount(2)+ CacheMissCount(3) + CacheMissCount(4)CacheMissCount(1) = 057CacheMissCount(2) N 3*DataTypeSizeBlockSizeN3CacheMissCount(3) dCacheMissCount(4) 2.5.N 2 DataTypeSizeBlockSizeВыводы ко второй главеВо второй главе представлена формула модели времени выполненияпрограмм для иерархии памяти.

Тем самым решается задача 1.В начале главы приводится абстрактная формула времени выполненияпрограмм. На основе данной формулы в параграфах 2.1-2.2 строятсяспециализированные формулы времени выполнения программ, реализующихалгоритм блочного умножения матриц, а также алгоритма Флойда-Уоршалла.

Впараграфе 2.3 описывается метод статического определения количества кешпромахов, возникших во время выполнения алгоритма.58Глава 3Умножение матриц рекордной производительностиНа протяжении последних 40 лет активно развиваются различные алгоритмыумножения матриц [17], [26].

Это, в частности, связано с тем обстоятельством, чток решению задачи умножения матриц сводятся алгоритмы решения других задач(например, один из алгоритмов решения СЛАУ [94], QR разложение матрицы[24], [25], LU-разложение матрицы [27]).На настоящее время различные алгоритмы умножения матриц реализованыво многих пакетах численных методов (Lapack, MKL, OpenBLAS, PLASMA ит.д.). Такое разнообразие библиотек обусловлено тем обстоятельством, чтосложно охватить все множество существующих вычислительных систем.Учитывая тот факт, что вычислительные системы с каждым годом стремительноусложняются, особый интерес представляют новые высокопроизводительныеалгоритмы, а также используемые ими методы оптимизации.Одним из главных факторов, сдерживающим достижение высокойпроизводительности, является скорость доступа к оперативной памяти, котораясущественно ниже скорости выполнения арифметических операций.

По этойпричине некоторые алгоритмы, которые использовались 20-30 лет назад [95], [96]сейчас практически не находят применения [76], [100]. В работе [71] показано,что алгоритм Штрассена, по сравнению со стандартным алгоритмом умноженияматриц, требует большее количество операций с памятью, что делает егонеприменимым на практике. Кроме того, очевидно, что он хуже по сравнению состандартным алгоритмом распараллеливается.В данной работе представлено описание нового алгоритма умноженияматриц, в основе которого лежит использование нестандартного размещения59матриц в оперативной памяти. В программной реализации задействованыразличные методы оптимизации для процессора, имеющего кеш-память иподдержку векторных вычислений. Исследуемая программа была протестированана процессоре с поддержкой 256-битных векторных регистров AVX.

Взаключительной части главы представлены результаты численных экспериментов,в которых проводилось сравнение реализованного алгоритма с функциейDGEMM из библиотек MKL, OpenBLAS, PLASMA.3.1.Высокопроизводительные пакеты программ линейной алгебрыСтатья [26] посвящена построению высокопроизводительного алгоритмаумножения матриц. В ней рассматривается семейство блочных алгоритмовумножения матриц и с помощью аналитических выкладок находит оптимальныйалгоритм, при котором время подкачки блоков сравнимо со временем умноженияблоков.

Исходная реализация алгоритма представлена в библиотеке GotoBLAS[86].На данный момент существует много высокопроизводительных библиотеклинейной алгебры. Среди них стоит выделить OpenBLAS [30], MKL [84],PLASMA [28], ATLAS [29].БиблиотекаOpenBLAS[30]являетсябиблиотекойсоткрытымпрограммным кодом, в основе которой лежит библиотека GotoBLAS. Онаоптимизирована под процессоры архитектуры Sandy Bridge, Haswell и Loongson.В частности, в функции умножения матриц данной библиотеки реализованаметодика упаковки регистров, позволяющая эффективно задействовать векторныевычисления.Библиотека PLASMA [28] является библиотекой решения задач линейнойалгебры,специальнооптимизированнойдляработынасовременныхмногоядерных процессорах. Библиотека PLASMA отличается от своих аналогов(OpenBLAS,MKL)тем,чтовсеалгоритмы,входящиевеесостав,оптимизированы для эффективного использования памяти.

Одной из причинвысокой производительности этих алгоритмов является то обстоятельство, что60они работают только с матрицами, размещенными в памяти блочным образом. Нанастоящее время поддерживаются алгоритмы решения нескольких наиболееважных задач линейной алгебры, таких как QR-разложение матрицы, LUразложение матрицы, умножения матриц. Следует однако отметить, что в этойбиблиотеке нет поддержки эффективной работы с разреженными матрицами(sparse matrices) и матрицами-полосами (band matrices).Библиотека ATLAS (Automatically Tuned Linear Algebra Software) являетсябиблиотекой с открытым исходным кодом.

Главной ее особенностью является то,что она не использует низкоуровневую оптимизацию. Таким образом, ее легчепереносить на новые платформы. В основе библиотеки ATLAS лежат алгоритмы,эффективно использующие память за счет блочных вычислений. Размеры блоковна конкретном компьютере подбираются эмпирически.3.2.Высокопроизводительный алгорит умножения матриц, использующийдвойное блочное размещение данныхБудем предполагать, что умножаются матрица A  M MxK и матрица B  M KxN ,а результат умножения сохраняется в матрицу C  M MxN . В основе программнореализованного автором алгоритма умножения лежит следующий блочныйалгоритм умножения матриц, представленный на листинге 1.Листинг 13.

Основа алгоритма умножения матриц рекордной производительностиvoid MatrixMult(double* A, double* B, double* C) {for (dk=0; dk<K; dk+=Dk)for (di=0; di<M; di+=Di)for (dj=0; dj<N; dj+=Dj)BlockAAddr = … // вычисление адреса блока матрицы ABlockBAddr = … // вычисление адреса блока матрицы BBlockCAddr = … // вычисление адреса блока матрицы CBlockMult(BlockAAddr, BlockBAddr, BlockCAddr);}61Внутри функции BlockMult производится умножение блоков матрицы A иматрицы B, результат сохраняется в блок матрицы C.

Матрицы при реализациитакогоалгоритмаумноженияразбиваютсянаблокисогласносхемепредставленной на Рисунке 7.Рис. 7: Схема разбиения матриц для высокоуровневой части алгоритмаумножения матрицЗначение Dk является одновременно шириной полосы матрицы A и высотойполосы матрицы B.

Для того, чтобы ускорить время подкачки данных изоперативной памяти в кеш-память, данные в матрицах изначально размещаютсяблочно. Матрица A размещается в оперативной памяти блоками размера Di .  Dk ,которые хранятся по столбцам (Рисунок 8).Рис. 8: Размещение блоков матрицы A в оперативной памяти.

Стрелкипоказывают, где блоки матрицы хранятся в оперативной памяти.62Соответственно, формула вычисления адреса блока матрицы A имеет вид:BlockAAddr(di, dk)  dk  M  di  DkМатрица B размещается в оперативной памяти блоками размера Dk .  D j ,которые хранятся по строкам (Рисунок 9).Рис. 9: Размещение блоков матрицы B в оперативной памяти. Стрелкипоказывают, где блоки матрицы хранятся в оперативной памяти.Соответственно, формула вычисления адреса блока матрицы B имеетследующий вид:BlockBAddr (dk, dj)  dk  N  dj  DkМатрица C разбивается на блоки размера Di .

 D j , которые хранятся построкам (Рисунок 10).63Рис. 10: Размещение блоков матрицы C в оперативной памяти. Стрелкипоказывают, где блоки матрицы хранятся в оперативной памяти.Соответственно, формула вычисления адреса блока матрицы имеет вид:BlockCAddr(di, dj)  di  N  dj  DiДля того, чтобы эффективно использовать кеш-память, размер блокаматрицы A должен принимать максимальное значение и помещаться целиком вкеш-память, в то время как D j  Min( Di .Dk ) . В реализованном алгоритме D j  8 .Алгоритм умножения блоков матрицы A и B является также блочным.Листинг 14. Умножение блоков в рассматриваемом алгоритме умножения матрицvoid BlockMult(double* BlockA,double* BlockB,double* BlockC){for (dk=0; dk<Dk; dk+=L)for (di=0; di<Di; di+=4)for (sk=0; sk<L; ++sk)for (si=0; si<4; ++si)for (sj=0; sj<8 ++sj)BlockC[di*8+sj*4+si] +=64BlockA[dk*Di+di*L+sk*4+si]*BlockB[dk*8+sk*8+sj];Умножаемые блоки при таком размещении разбиваются на блоки меньшегоразмера согласно схеме, представленной на Рисунке 11.Рис.