Разработка модуля расчета моделируемой схемы с разреженными матрицами

Таким образом, вернемся к упомянутому ранее разреженному строчному формату хранения и наглядно продемонстрируем преимущества такого подхода, приведя сопоставление самых распространенных форматов хранения матриц разреженного типа (двумерный массив, координатный формат (COO), строчный формат (CRS)). Результаты этого сравнения указаны в табл. 2.



Таблица 2 -- Сравнение избыточности форматов хранения разреженных матриц

Размерность разреженной матрицы	Количество памяти (МБ) выделяемой под хранение разреженной матрицы целочисленных значений для каждого формата
	Двумерный массив	COO	CRS
1000 Ч 1000	3.814	1.831	1.224
5000 Ч 5000	95.367	45.776	30.536
12000 Ч 12000	549.316	263.671	175.827
20000 Ч 20000	1525.878	732.421	388.357

График зависимости количества выделенной памяти от размерности матрицы и выбранного формата хранения представлен на рис. 15.

Рисунок 15 -- График зависимости количества выделенной памяти от размерности матрицы и выбранного формата хранения

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

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

Ключевой модификацией метода, реализуемого в данной работе, является использование CRS формата. За счет применения CRS формата для хранения промежуточных итогов вычисления удалось значительно повысить уровень оптимизации производительности, что будет продемонстрировано в главе 4. Кроме того, благодаря этому происходит существенное снижение времени выполнения разработанного приложения.



4. Последовательный алгоритм умножения разреженных матриц

Помимо способа хранения матриц, следует обсудить подходы к расчету возникающих разреженных систем. Существует ряд подходов к обработке и вычислению разреженных матриц, которые способны тем или иным образом увеличить общую скорость алгоритма.

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

Таблица 3 Время работы алгоритма умножения матриц при разных размерах блока

Размер блока	16	32	64	128
N=1024	14.3 с	9.318 с	6.058 с	6.606 с

Рисунок 16 -- Диаграмма времени работы алгоритма умножения матриц при разных размерах блока

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

Говоря о действиях с матрицами, в частности, обсуждая проблемы матричного умножении, следует упомянуть асимптотическую сложность алгоритмов. Для оценки сложности алгоритмов была введена матричная экспонента щ -- число с таким минимальным значением, что при любом е > 0 и любом достаточно большом количестве n число операций арифметического характера, требуемых для того, чтобы вычислить произведение двух матриц, имеющих размер nЧn, будет не больше nщ + е.

Легко продемонстрировать, что стандартный алгоритм, используемый для перемножения матриц (предусматривающий сложение результатов произведений расположенных в одной из строк матрицы А элементов на аналогичные, помещенные в столбце матрицы В), требует совершения n3 умножений и n2 ?(n-1) сложений. На протяжении длительного периода он считался наиболее быстрым, однако в конце 60-х годов прошлого века появилась новая методика с щ < 2,808. Немецкий математик Штрассен предложил алгоритм с меньшим количеством совершаемых умножений (7 вместо 8) в матричной алгебре М2Ч2. Однако за счет рекурсии он успешно может быть использован для матриц более высоких порядков, при этом шаг рекурсии сопровождается повторением одинаковых семи перемножений. Минусом способа Штрассена является его относительная сложность и необходимость использования большого объема памяти при программировании.

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

Рисунок 17 -- График изменения щ

Пара американских ученых Копперсмит и Виноград сумели довести показатель до щ < 2.376, и после этого, то есть с начала 90-х годов прошлого века, оценка матричной экспоненты снижаться перестала. Теоретически все это время именно данный алгоритм был самым быстрым, но практически он не применялся, поскольку его эффективность заметна лишь при условии, что матричные размеры сверх велики.

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

Целевого эффекта вследствие трехкратного исполнения алгоритма получено не было, а оценивать щ при новых условиях стало существенно труднее. Вследствие этого рядом ученых математиков был сделан вывод о том, что достигнутое на заре 90-х значение является итоговым.

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

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

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

Существует две разновидности метода последовательной реализации. Первый способ основан на двойном транспонировании матрицы и выделении численного и символического этапа. Посредством данного алгоритмического подхода осуществляется матричное преобразование не только из строчного вида в столбцовый, но и из столбцового в строчный. Так как матрица показанная в строчном варианте является матрицей в столбцом варианте по прошествии фазы транспонирования, то по сути можно утверждать, что данный алгоритмический подход реализует действие транспонирования. Интересной особенностью подхода является тот факт, что итоговое представление выражено в упорядоченном виде, так как столбцовые индексы элементов по строкам находятся в порядке возрастания. В этом смысле, если применить к изначально неупорядоченной матрице два раза операцию транспонирования озвученного алгоритма, то результатом будет будет та же матрица, но в упорядоченном варианте. Для матрицы симметричного типа, у которой полное представление, достаточно лишь одного транспонирования для получения упорядоченного варианта. Однако следует учитывать, что в случае с симметричными матрицами чаще прибегают не к полному представлению, а к верхнему. Задача такого представления разреженной матрицы (А) может быть рассмотрена с позиции, применяемой для матриц общего вида. При этом заданное представление будет просто восприниматься в качестве полной проекции другой матрицы (В). Области верхних треугольников двух данных объектов совпадают. При этом нижний треугольник в B, наряду с главной диагональю, содержит нулевые значения. При таких условиях для упорядочения В достаточно двойного транспонирования, благодаря которому в результате будет получен упорядоченный вариант верхнего треугольника А.

Следовательно, для того, чтобы упорядочить матрицу, которая задана в верхнем представлении, необходимо провести двойное транспонирование, направленное лишь на половину матрицы. В связи с этим, когда матрица представлена в полном варианте, можно осуществить только одно транспонирование. Немного изменив представленный алгоритмический подход, ему возможно найти применение в вопросе перестановки столбцов/строк матрицы разреженного типа [23]. Интересная возможность возникает в случае, когда какая-либо операция над матрицей осуществляется с помощью алгоритма, который включает в себя численную и символическую фазы. Еще лучше, когда результат является упорядоченным. Так как матричный портрет, который определен символической фазой, не будет изменен численной фазой, то наиболее разумным кажется упорядочивание лишь портрета матрицы до этапа расчетов.

Действия запланировано осуществлять в определенной последовательности:

Задействуя символическую фазу озвученного подхода, требуется найти портрет итоговой матрицы

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

Конечным шагом будет расчет значений элементов посредством численной фазы

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

Другим интересным моментом можно назвать исключение Гаусса. Его возможно сделать лишь по прошествии фазы упорядочивания итогового матричного портрета. Опираясь на тот факт, что операции осуществляются с матрицей разреженного типа, следует применять исключение Гаусса по строкам. Улучшение показателей эффективности можно получить с помощью задействования матричных элементов в последовательности, в которой они находятся в матрице, либо другими словами в естественной последовательности. Исключение Гаусса, скажем, у матрицы А, осуществляется в n этапов по строкам. В исходной точке этапа k есть матрица A(k), в которой слева от диагонали в строках от первой до (k-1)-ой имеются элементы со значением 0. Также в матрице находятся значения 1 в первых (k-1) диагональных положениях. Далее можно увидеть внешний вид A(k) матрицы при k = 4, n = 6:

Рисунок 18 -- Структура матрицы А(k)

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

Нужно подчеркнуть три этапа озвученного алгоритмического подхода:

деление каждого элемента матрицы, который расположен в 1-ой строке на находящийся в положении (1. 1) элемент диагонали.

исключение значения, которое находится в положение (2. 1). Чтобы это сделать производится вычитание из строки под номером 2 произведения первой строки и данного элемента. Далее строка 2 проходит через операцию нормировки, посредством деления всех ее значений на элемент диагонали из положения (2. 2).

Заключительным этапом является исключение значения, которое находится в положении (3. 1).Чтобы осуществить данное действие нужно вычесть из строки под номером 3 произведение строки 1 и требуемого скаляра. После этого элемент, который находится в положение (3. 2) тоже исключается, для чего осуществляется порядок действий, идентичный озвученному, с той лишь разницей, что вместо строки под номером 1 происходит подстановка строки 2 и тд.

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

Те значения, которые находятся справа в строке 3, могут быть подсчитаны посредством отнятия из 3-ей строки первых двух строк, пониженных на соответствующие множители данной строки 3, после чего последняя нормируется. Различие двух алгоритмов заключается только в порядке, в соответствии с которым делается исключение. Но при этом должны также приниматься во внимание вопросы скорости.

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

Пусть массивы под названием IA, JA, AN являют собой некое представление разреженной матрицы размера NхM. Необходимо добиться создания массивов IAT, JAT, ANT, которые хранят разреженное представление после операции транспонирования. Желая узнать, что за элементы одновременно имеются в IA, JA, AN и также располагаются в данном столбце, мы приходим к неэффективному проходу по JA и IA. Чтобы улучшить ситуацию, требуется использовать M целых перечней. Они в самом начале будут пустыми. Также необходимо отметить для всех из них 1-ое положение, которое свободно. В начале, очевидно, будут использоваться 1-ые положения. Помимо этого нужно сделать инициализацию M списков вещественного типа. Далее осуществим просмотр каждого значения матрицы, которое имеет значение, отличное от нуля. Для всех тех, которые располагаются в IA, сделаем добавление текущего элемента в J-ый перечень, а значение данного элемента в J-ый вещественный перечень. После этого осуществляем операцию инкрементирования значения нужного указателя. Покажем демонстрацию озвученных действий. Осуществим анализ представленной разреженной матрицы:

Рисунок 19 -- Схематичный вид матрицы разреженного типа

Подход RR (C)U для нее представлен в следующем виде:

Здесь не отображен AN. В связи с тем, что матрица насчитывает шесть столбиков, требуется инициализировать 6 списков. После чего нужно пройти строку под номером 1, имеющую индексы столбцов 5, 6 и 3; в связи с этим нужно интегрировать в озвученные списки единицу:



После этого располагаем двойки в списках под номером 4 и 1:

Далее, по прошествии обработки все 5-и строк, обсуждаемые списки будут представлены следующим образом:

или по другому:

Представленные рассуждения, вплоть до массивов ANT и IAT, являются подходом CR (C) O заданной матрицы, или схемой RR (C) O, но уже матрицы, которая прошла операцию транспонирования. IAT находится посредством элементарного вычисления значений в каждом из списков. ANT при необходимости можно найти, интегрируя значение всех отличных от нуля элементов, находящихся в AN, в список вещественного типа в тот момент, когда необходимое число целого типа вставлено в целый список. Касаясь практической стороны вопроса, списки могут быть созданы прямо в массивах JAT и ANT, а указатели, которые хранят 1-ое незанятое положение каждого из списков - в IAT. Из этого следует, что не нужно использовать добавочную память, а перемещение данных, которые находятся в памяти компьютера, минимизировано.

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

Рисунок 21 -- Вариант портрета разреженной матрицы

Рисунок 22 -- Вариант портрета разреженной матрицы

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

Стоит обратить внимание на тот факт, что обычно строки в формате CRS рассматриваются по порядку, однако значения в строках могут являться как упорядоченными по номеру столбца, так и неупорядоченными. То, как элементы расположены, легко понять, в связи с тем, что массив под названием Col содержит в себе номера столбцов. У обоих подходов есть своя собственная специфика по отношению к последующей фазе обработки. Оба варианта, как упорядоченный, так и неупорядоченный, могут быть использованы.

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

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



Рисунок 23 -- Визуализация умножения матриц

Исходя из определения, получается, что принадлежащий к строке i и столбцу j элемент матрицы C представляет собой скалярное произведение i-ой матричной (А) строки и j-го матричного (В) столбца. Рассмотрим ключевые особенности, связанные с разреженным представлением.

Первая из них состоит в том, что задействуемая структура данных, в основе которой лежит формат CRS, предусматривает хранение исключительно тех элементов, которые имеют ненулевое значение. Это делает программирование вычисления скалярного произведения более сложным процессом, однако, вместе с тем, число требуемых арифметических операция снижается. Рассчитывая скалярные произведения, не требуется осуществлять перемножение нулей и накопление полученного нуля в частичную сумму. Это отражается на времени, затрачиваемом на подсчет, в сторону его уменьшения.

Допустим, в первом и втором векторах присутствует по 1% элементов с ненулевым значением, но лишь каждый десятый из них занимает позицию, которой соответствует. При таких исходных данных расчет с применением сведений относительно структуры векторов может быть произведен с задействованием умножений и сложений, число которых 1000-кратно меньше. С учетом того, что подобных комбинаций в квадрате N-ное количество, объем вычислений кардинально сокращается. Но на практике все оказывается несколько сложнее. Чтобы структура векторов была должным образом учтена, также необходимо машинное время. Приходится сопоставлять номера элементов с ненулевым значением, чтобы обнаруживать пару значений, которые должны быть перемножены и накоплены в частичную сумму. Вследствие этого алгоритм начинает ветвиться.

Вторая особенность заключается в том, что необходимо в матрицах А и В научиться определять векторы. Исходя из определения, имеются в виду строки матрицы А и столбцы матрицы В. Легко произвести выделение матричной строки в формате CRS: для i-ой строки уже известны исходные данные в виде первого элемента (RowIndex[i]) и последнего (RowIndex[i+1]-1), благодаря чему можно узнать значения элементов и номера, присвоенные столбцам, которые хранятся в массивах соответственно Value и Col. В итоге на проход по строке расходуется время, пропорциональное количеству элементов с ненулевым значением, которые внесены в указанную строку, на проход по всей совокупности строк - время, которое пропорционально числу элементов с ненулевым значением NZ.

Сложность состоит в том, чтобы выделить столбец. Для нахождения элементов столбца j требуется выполнить просмотр массива Col (с приблизительным количеством операций, равных NZ) и обозначить все множество элементов, которые имеют число j в соответствующей ячейке массива Col. При необходимости выделить ту же операцию для всех столбцов требуется ориентировочно NZ*N операций. Данный путь является неэффективным (N - порядок матрицы).

Выйти из ситуации можно за счет последовательного транспонирования матрицы В. При достаточно быстром исполнении это даст возможность работать со столбцами стартовой матрицы В с высокой степенью эффективности.

Третья особенность заключается в том, что требуется умение записывать элементы, которые были подсчитаны, в матрицу С. При этом важно не допускать перепаковки, помня о том, что формат хранения матрицы С - CRS. Поэтому следует позаботиться о том, чтобы матрица С пополнялась элементами, имеющими ненулевое значение, в строгой последовательности, построчно, начиная с верхней и продолжая вниз, двигаясь при этом слева направо. Отсюда следует, что нередко бумажный метод вычисления, предполагающий написание первого столбца матрицы В над матрицей А с последующим перемножением всех строк, не является приемлемым. Поступать следует как раз наоборот, беря первую строку матрицы А и умножая ее на каждый столбец матрицы В поочередно. При таком подходе матрица С заполняется последовательно, что дает возможность вносить элементы в массивы Value и Col, подвергать форматированию массив RowIndex непосредственно в процессе.

Итак, следовать необходимо такому алгоритму:

Создать алгоритм транспонирования разреженной матрицы и применить данный подход к матрице В.

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

Перемножить последовательно все строки из матрицы А на все строки матрицы BT, внося результаты в C и формируя тем самым ее структуру. При умножении строк реализовать сопоставление с целью выделения пар ненулевых элементов.

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

В чем суть операции транспонирования? Если создать нулевую матрицу, а далее добавлять туда по одному элементу, выбирая их из CRS-структуры исходной матрицы, придется столкнуться с необходимостью большого количества перепаковок. Чтобы этого избежать, нужно формировать транспонированную матрицу построчно. Для этого можно брать столбцы исходной матрицы и создавать из них строки результирующей матрицы. Однако выделить из CRS-матрицы столбец i не так просто. Необходимо другое решение.

Обнулим массив RowIndex матрицы АТ. Просмотрим массив Col матрицы A и сосчитаем количество элементов в каждом столбце матрицы A (строке матрицы АТ), сохраняя результаты в массиве RowIndex матрицы АТ. Пусть при этом AT.RowIndex[j] хранит количество элементов в столбце j-1 матрицы А (j меняется от 1 до N включительно).

Рисунок 24 -- 1 этап транспонирования

Подсчитаем индексы начала каждой строки в матрице АТ так, что элемент AT.RowIndex[i] хранит индекс начала (i-1)-ой строки.

Рисунок 25 -- 2 этап транспонирования

Правильно расположим в структуре результирующей матрицы значения из исходной матрицы. Сами значения и их номера строк (теперь столбцов) известны, необходимо «лишь» поместить их в корректную позицию. Для этого будем использовать массив AT.RowIndex. Заметим, что взяв очередной элемент матрицы A, имеющий значение V, координаты (i, j), мы должны добавить его в j-ю строку матрицы AT. Зная, что в настоящий момент j-я строка матрицы AT начинается с элемента AT.RowIndex[j+1], будем добавлять V и i в массивы AT.Value и AT.Col соответственно по адресу AT.RowIndex[j+1], после чего увеличим AT.RowIndex[j+1] на единицу. В конце элемент AT.RowIndex[j+1] будет хранить индекс начала (j+1)-ой строки, что и требуется.

Рисунок 26 -- 3 этап транспонирования

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

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

Также стоит обсудить в целом алгоритм умножения разреженных матриц [24], который реализован в рамках алгоритма метода последовательной реализации.

В начале происходит передача в функцию требуемых массивов данных, а именно: массивов обычных и транспонированных матриц формата CRS в виде векторов, также передаются массивы выходных значений для результирующей матрицы. После этого происходит выделении памяти в зависимости от размеров матрицы, после чего создается временной массив temp, который нацелен на хранение промежуточных результатов расчета. Далее следует непосредственно алгоритм работы. Внешний цикл ходит по строкам первой матрицы, внутренние циклы в свою очередь обрабатывают элементы строк, а также строки второй матрицы с их элементами. Первый внутренний цикл “идет” по элементам текущей строки первой матрицы и запоминает столбец, а также индекс начала строки данного элемента. Второй внутренний цикл “проходит” по строкам второй матрицы, а третий также запоминает столбец текущего элемента, а потом проверяет, что из первой матрицы элемент для того же столбца существует (не равен 0), чтобы не осуществлять лишние действия и впоследствии производить умножение и сложение элементов. Как только данная последовательность действий была проделана для всех элементов текущей строки второй матрицы, происходит заполнение выходных массивов CRS. После прохождения по всем строкам второй матрицы осуществляется изменение массива индексов начала строки в соответствии с требуемым количеством элементов. Озвученный алгоритм повторяется для всех строк первой матрицы. 

Рисунок 27 - Часть алгоритма умножения



5. Параллельный алгоритм умножения разреженных матриц

Традиционная ЭВМ была последовательной, то есть в определенный момент времени только над одним операндом выполнялась только одна операция.

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

Параллельные вычисления разделяются на две парадигмы: message passing или параллелизм задач и data parallel или параллелизм данных. За фундамент этих подходов приняли идею распределения вычислительной работы по процессорам доступным пользователю. Параллелизм данных характеризуется:

одна программа управляет обработкой всех данных;

пространство имен является глобальным (детали скрыты от программиста);

на параллельных процессах происходит слабая оптимизация вычислений;

на всех доступных данной программе процессорах параллельные операции выполняются одновременно.

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

Параллелизм задач характеризуется:

минимизация обмена данными между задачами;

повышение трудоемкости программы и ее отладки;

повышенная опасность возникновения тупиковых ситуаций;

ответственность за разгрузку процессоров лежит на программисте.

Перед тем как перейти к классам ЭВМ давайте уточним способы решения задач поставленных перед ЭВМ методами параллельных вычислений. Можно выделить два основных метода: конвейер и параллелизм.

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

система работает над выполнением повторяющейся операции;

выполняемую операцию можно разделить на самостоятельные части;

трудоемкость подоперации примерно одинакова.

Рисунок 28 - ЗУ - Запоминающее устройство; Пр - процессор

В отправляемые на конвейер данные должны быть независимы. Если же данное условия не выполняются и существуют операнды, зависящие от предыдущих, то в конвейере образовываются так называемые “пузыри”. Это еще одна проблема в работе конвейерных систем.

Параллельное решение задачи предусматривает распределение выполняемых команд и операндов на несколько процессоров [26].

Рисунок 29 - ЗУ - Запоминающее устройство; Пр -процессор

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

5.1 Классификация Флинна. Архитектура вычислительного оборудования

Существуют четыре основных класса вычислительного оборудования по потокам данных и команд:

SISD;

MISD;

SIMD;

MIMD.

В первом случае (SISD - single instruction single data) представлена последовательная вычислительная машина. Одна команда выполняется один раз над одним элементом данных. Сюда относятся все машины с одним процессором и одним комплектом памяти.

Далее идет MISD (multiple instruction single data) - множественный поток данных выполняется над одним элементом данных. Архитектура с единообразной памятью. Сюда относятся многоядерные машины, машины с конвейерными архитектурами процессоров, симметричные многопроцессорные, многопроцессорные графические адаптеры.

В SIMD (single date multiple instruction) один поток команд выполняется над множеством элементов данных. Этот тип характерен для векторных машин и для машин с параллельными данными. В структуре можно выделить множество вычисляющих устройств. Сопоставление команд происходит векторным компилятором. Данный класс приспособлен к однообразным вычислениям.

В MIMD (multiple data multiple instruction) множество команд выполняется множественным потокам данных. Для того, чтобы задача была выполнена, необходимо обеспечить синхронизацию и обмен информацией между узлами, входящими в вычислительную систему. К данной категории причисляются кластеры и грид-сети, пребывающие в режиме обработки.

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

Посредством разделяемой памяти

На каждом процессоре, входящем в многопроцессорную систему, происходит запуск исполнительного потока, относящегося к одному процессу. Обмен данными между потоками осуществляется благодаря наличию общего участка памяти, выделенного под определенный процесс. Число потоков равно числу используемых в системе процессоров. Создание потока осуществляется посредством языков программирования или через библиотеки с применением компиляционных средств. Традиционно такое параллельное программирование нуждается в наличии некой формы захвата управления, чтобы налаживать межпоточную координацию.

Посредством передачи сообщений

Производится запуск однопоточного процесса на всех процессорах, составляющих комплексную многопроцессорную систему. Осуществляется обмен информацией между разрозненными процессами, за течение которых отвечают разные процессоры, посредством специальных сообщений. Создание процессов происходит в явной форме, для чего вызываются необходимые функции операционной системы, для передачи же сообщений задействуется библиотека или средства языка программирования. Предусматривается возможность отправки сообщений в асинхронном режиме или в режиме «рандеву», при котором сторона, отправившая сообщение, оказывается заблокированной до того момента, пока адресат его не получит. Также существует разделение асинхронной передачи на надежную и ненадежную.

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

Характеристики производительности существенно разнятся, когда речь идет о параллелизме с памятью распределенного типа и с передачей сообщений. Преимущественно системы с передачей сообщений отличаются более низкими накладными расходами памяти на процесс и времени на смену задач, хотя при этом больше ресурсов уходит собственно на передачу сообщений по сравнению с тем, что требуется для вызова процедуры. Нередко другие факторы накладывают серьезный отпечаток на производительность, перекрывая приведенные различия.

Гибридный способ

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

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

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

5.2 Сравнение технологий распараллеливания

Чтобы выполнить поставленную задачи необходимо определить наиболее эффективный алгоритм. Для сравнения были выбраны четыре технологии распараллеливания вычислений: автоматическое распараллеливание, Open MP, MPI, CUDA. Для этой задачи покажем теоретические и расчётные ускорения, получаемые технологиями MPI, OpenMP и автоматическим распараллеливанием программ компилятором GCC с использованием инструкций векторных вычислений SSE2 и AVX2.

В данной работе в качестве задачи, решаемой с использованием различных технологий параллельных вычислений, будет использоваться умножение квадратных матриц методом Фокса. Матрицы - A , B и C=AB . При таком методе решения используется блочное представление матриц так, что последнее равенство можно представить в следующем виде:



Количество процессоров одновременно решающих задачу равно , где . За нахождение блока матрицы ответственен процессор :

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

Заметим, что обмен данными будет широко применяться при использовании первой технологии - MPI. При решении задачи технологией OpenMP явного обмена между процессорами нет ввиду использования общей памяти. Наконец, при автоматическом распараллеливании вся задача будет решаться одним процессором так, что размер исходной матрицы будет равен размеру единственного выделяемого из такой матрицы блока, получая ускорение за счёт использования процессом векторных инструкций, обрабатывающих за команду набор данных.

Message Passing Interface (MPI)

Параллельной называют программу в рамках MPI, которая представляет собой большое количество запущенных и одновременно выполняемых процессов. Это может происходить на обособленных процессорах либо на одном, но при условии разбивки из реализации по времени. Параллельная программа в редких случаях может реализовываться на одном процессоре. Обычно к такому способу прибегают, когда требуется выполнить стартовую проверку, насколько правильно функционирует параллельная программа [27].

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

При старте параллельной программы, запускаемой за счет средств исполнения MPI, определяется, сколько процессоров будет использоваться, и сколько процессов они будут обеспечивать. Первоначальное количество не подлежит изменению в процессе вычислений без специальных средств, направленных на порождение и динамическое управление ими. Каждый процесс и программа имеют последовательную нумерацию, при этом номер, присваиваемый процессу, называется его рангом.

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

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

Коммуникатор - еще одно значимое для MPI понятие, которое описывает комплекс процессов. Данный термин употребляется по отношению к специально формируемому служебному объекту, объединяющему в себе совокупность процессов и множество дополнительных параметров, которые необходимы для осуществления операций передачи информации.

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

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

В MPI присутствует широкий набор данных, относящихся к базовым типам, которые в значительной степени аналогичны тем, которыми располагают алгоритмические языки C и Fortran. Также важно отметить, что MPI содержит возможности, позволяющие формировать новые производные типы данных, которые позволяют четче и короче описывать контент пересылаемых сообщений.

Таким образом, MPI можно назвать библиотекой функций, программным интерфейсом (API), ориентированным на системы, имеющие распределенную память для трансфера информации, и позволяющим производить обмен сообщениями между нацеленными на реализацию одной задачи процессами.

Данная библиотека располагает чрезвычайно широким функционалом, обеспечивающим, помимо прочего:

реализацию операций парного типа;

инициализацию и закрытие MPI-процессов;

реализацию операций коллективного характера;

работу со структурными данными;

работу с группами процессов и коммуникаторами;

формирование топологии процессов.

Рисунок 29 - Схематичное изображение процесса выполнения программы с использованием MPI

Теоретическое время выполнения программы на p процессорах:

,

, (7),

где n - размер матриц, - время инициализации, - время выполнения одной арифметической операции, - время пересылки.

Теоретическое ускорение, получаемое при p процессорах:

, (8),

Значения параметров - , с, с, c.

При вычислении расчётного ускорения использовалось среднее значение трёх измерений. На следующем графике показано сравнение теоретического и расчётного ускорений:

Рисунок 30 - Сравнение теоретичесого и расчётного ускорения MPI

Рисунок 31 - Реализация умножения с использованием MPI

Open Multi-Processing (Open MP)

OpenMP (Open Multi-Processing) - это открытый стандарт для языков C, C++ и Fortran, позволяющий распараллеливать программы, описывающий компиляторный комплекс директив и переменных окружения, участвующих в программировании приложений многопоточного вида посредством мультипроцессорных систем, укомплектованных совместной памятью.

Для распараллеливания в OpenMP используется прием интеграции в программный текст спецдиректив, вызова функций вспомогательной категории и применения переменных окружения. Для OpenMP характерна SPMD-модель программирования параллельного типа, предусматривающая использования одинакового кода для всех параллельных нитей [28].

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

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

Когда все параллельные потоки выполнены, программа снова переходит под управление главного потока. На данном этапе важно правильно организовать трансфер данных от параллельных потоков основному. В данном контексте на первый план выходит синхронизация завершения работы параллельных потоков, так как нередко даже те процессы, трудоемкость которых сопоставима, могут иметь совершенно различное время выполнения [29-30]. В OpenMP процесс обмена данными осуществляется посредством общих переменных. В связи с этим возникает потребность разграничивать одновременный доступ разных нитей к пребывающим в общем пользовании данным. Достигается это за счет развитых синхронизационных средств. Важно помнить, что чрезмерное увлечение синхронизациями оборачивается ощутимым замедлением программы.

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

Рисунок 32 - Схематичное изображение процесса выполнения программы с использованием OpenMP

Теоретическое время выполнения программы на процессорах:

, (9)

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

Теоретическое ускорение, получаемое при процессорах:

, (10)

Функция - аппроксимация роста доступа к памяти при конкуренции процессоров, при доступ единственного процессора не изменится (), но при увеличении числа процессоров - возрастёт (). График функции :



Рисунок 33 - График аппроксимации роста доступа к памяти

Ввиду того, что , размер элемента матрицы (int) - 4 байта, размер одной матрицы равен Мб. Доступ к элементам матрицы - не последовательный, в связи с чем, матрица не помещается целиком ни в один кэш, поэтому время доступа к оперативной памяти исчисляется от 100 до 1000 тактов, поэтому положим равным .

Значения параметров - , с, с, , .

При вычислении расчётного ускорения использовалось среднее значение трёх измерений. На следующем графике показано сравнение теоретического и расчётного ускорений:

Рисунок 34 - Сравнение теоретичесого и расчётного ускорения OpenMP



Рисунок 35 - Реализация умножения с использованием OpenMP

Автоматические распараллеливание

Автоматическое распараллеливание -- способ автоматического преобразования программы, которое осуществляет компилятор для параллелизации вычислений.

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

Данный вид распараллеливания программ нельзя назвать эффективным, так как здесь приемлимые показатели ускорения могут быть достигнуты только для простых циклов for, где отсутствуют зависимости по данными между итерациями и при этом число итераций не может измениться после начала цикла. При условии, что два указанных условия в некотором for-цикле выполняются, но он имеет сложную неочевидную структуру, то он не будет распараллеливаться.

Виды автоматического распараллеливания:

Полностью автоматический: компилятор выполняет распараллеливание самостоятельно, без вмешательства программиста.

Полуавтоматический: специфические аспекты организации распараллеливания задаются программистом с помощью специальных ключей. Автоматическое распараллеливание имеет следующие недостатки:

Есть вероятность возникновения ошибочное преобразование логической составляющей программы

Есть вероятность увеличения времени выполнения программы

Нет гибкости ручного распараллеливания

Только циклы распараллеливаются эффективно

Если алгоритм сложный, то распараллеливание может не работать

В качестве будет рассматриваться выполнение программы без использования автоматического распараллеливания, при - распараллеливание с инструкциями SSE2, обрабатывающими за такт 128 бит (128/32 = 4, 32 бита - размер элемента матрицы), при - инструкциями AVX2, использующими 256-битные регистры (256/32=8).

Теоретическое время выполнения программы на процессорах:



, (11)

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

, (12)

Ввиду организованного последовательного доступа к памяти, для обеих матриц и будет использоваться кэш процессора 1-го уровня (L1), время доступа к которому равно порядку нескольких тактов, поэтому .

Значения параметров - , с, с, , .

При вычислении расчётного ускорения использовалось среднее значение трёх измерений. На следующем графике показано сравнение теоретического и расчётного ускорений:

...