CC BY
70
ISSN 1992-6502 (Print)_
2015. Т. 19, № 3 (69). С. 178-189
Ъюьшм QjrAQllQj
ISSN 2225-2789 (Online) http://journal.ugatu.ac.ru
УДК 519.612
Двухуровневый параллельный алгоритм выполнения
численной фазы разложения Холецкого для разреженных матриц
12 3
С. А. Лебедев , И. Б. Мееров , Е. А. Козинов , Д. Р. Ахмеджанов4, А. Ю. Пирова5, А. В. Сысоев6
1sergey.a.lebedev@gmail.com, 2meerov@vmk.unn.ru, 3evgeniy.kozinov@gmail.com, 4dmitry.akhmedzhanov@yandex.ru, 5anna.yu.malova@gmail.com, 6sysoyev@vmk.unn.ru
16Нижегородский государственный университет им. Н. И. Лобачевского
Поступила в редакцию 10 августа 2015 г.
Аннотация. Рассматривается задача распараллеливания численной фазы разложения Холецкого для разреженных симметричных положительно определенных матриц. Предложена новая схема распараллеливания мультифронтального метода для систем с общей памятью. Данная схема основана на сочетании двух подходов к организации параллелизма на разных уровнях дерева исключения. В нижней части дерева выполняется параллельная обработка узлов, хранящихся в приоритетной очереди. На верхних уровнях дерева узлы обсчитываются последовательно с использованием многопоточного BLAS. Приведены результаты вычислительных экспериментов на матрицах из коллекции университета Флориды. Показана сопоставимость выполненной реализации с решателями MUMPS и MKL PARDISO.
Ключевые слова: алгебра разреженных матриц; разложение Холецкого; численная фаза; мультифрон-тальный метод; высокопроизводительные вычисления; динамическая схема распараллеливания; логические задачи.
ВВЕДЕНИЕ
Системы линейных алгебраических уравнений (СЛАУ) с разреженной симметричной положительно определенной матрицей возникают при моделировании процессов во многих предметных областях. Огромная размерность таких систем (в современных приложениях - 107 и выше) приводит к большим затратам памяти и процессорного времени, что без сомнения позволяет отнести их решение к области применения высокопроизводительных вычислений. Прямые методы, основанные на факторизации матрицы системы, активно применяются для решения больших разреженных СЛАУ. В настоящее время в этой области разработан целый ряд параллельных алгоритмов для систем с различной архитектурой и соответствующие про-
Работа рекомендована программным комитетом международной научной конференции «Суперкомпьютерные дни в России». Работа выполнена при частичной поддержке гранта РФФИ №14-01-31455, гранта МОН РФ (соглашение от 27 августа 2013 г. № 02.В.49.21.0003 между МОН РФ и ННГУ).
граммные пакеты [1], среди которых широко распространены MKL PARDISO, MUMPS, SuperLU, CHOLMOD и другие.
С 2011 г. на факультете ВМК ННГУ разрабатывается прямой решатель разреженных СЛАУ с симметричной положительно определенной матрицей, основанный на методе Холецкого [2]. В рамках данного решателя изучаются вопросы оптимизации соответствующих алгоритмов под современные многоядерные архитектуры. В данной работе рассматривается задача распараллеливания наиболее затратной по времени и памяти численной фазы разложения Холецкого. Несмотря на наличие множества алгоритмов и их реализаций, вопрос о развитии существующих методов построения масштабируемых алгоритмов и программных средств, ориентированных на системы с общей памятью, не потерял свою актуальность в связи с постоянным развитием многоядерных архитектур. Ранее в работе [3] мы предложили способ распараллеливания мультифронтального метода, основанный на динамической схеме балансировки нагрузки. В статье предлагается модификация, состоящая в сочетании возможностей базовой схемы с использованием параллельного
BLAS для решения «тяжеловесных» задач на верхних уровнях дерева исключения. Будет показано, что применение данной комбинированной схемы позволяет улучшить масштабируемость программной реализации.
ПОСТАНОВКА ЗАДАЧИ И МЕТОД РЕШЕНИЯ
Дана система линейных уравнений Ах = Ь, где А - разреженная симметричная положительно определенная матрица; b - плотный вектор; х- вектор неизвестных. Необходимо найти решение системы х.
Принцип работы прямых методов основан на факторизации матрицы системы с последующим решением треугольных систем. Для симметричных положительно определенных матриц факторизация выполняется методом Хо-лецкого. Для этого в большинстве случаев используется двухфазный подход: вначале находится портрет фактора, т.е. расположение ненулевых элементов (символьное разложение), затем полученный портрет заполняется значениями (численное разложение). Необходимо отметить, что символьная фаза выполняется гораздо быстрее численной, поэтому множество усилий исследователей направлено на оптимизацию и распараллеливание численной фазы.
Существует несколько методов выполнения численной фазы разложения Холецкого. Среди них можно выделить три наиболее широко используемых на практике: ориентированный влево (left-looking) [4], ориентированный вправо (right-looking) [4] и мультифронтальный (multifrontal) [5]. Основное отличие между ними заключается в способе формирования результирующей матрицы L, а также в способе хранения и размещении в памяти промежуточных результатов. Перечисленные методы показывают схожую производительность на различных тестовых наборах матриц, но с точки зрения многих исследователей мультифронтальный метод является наиболее перспективным для распараллеливания [6, 7]. По этой причине мультифрон-тальный метод используется в качестве базового в данной работе.
КРАТКОЕ ОПИСАНИЕ МУЛЬТИФРОНТАЛЬНОГО МЕТОДА
Впервые мультифронтальный метод был представлен Даффом и Рейдом [8] в 1983 г., а затем развит Лю [9], Амстоем и Даффом [10]. К числу основных достоинств мультифронталь-ного метода относят эффективное использова-
ние кэш-памяти всех уровней. При реализации с использованием подходящих структур данных основными становятся операции с плотными подматрицами, для выполнения которых может быть использован BLAS 3. Данный метод используется в одном из широко распространенных решателей c открытым исходным кодом MUMPS. К числу недостатков мультифронталь-ного метода можно отнести высокие затраты памяти для представления промежуточных результатов и большое число операций с плавающей точкой, особенно для задач, полученных путем дискретизации трехмерного пространства [11]. Приведем краткое описание метода.
Численная фаза разложения Холецкого применяется для заполнения уже сформированного шаблона фактора матрицы численными значениями. Для этого в мультифронтальном методе процесс факторизации разбивается на факторизацию небольших плотных подматриц, называемых фронтальными (frontal matrix). При этом порядок получения столбцов фактора определяется графом задач, который в случае симметричной положительно определенной матрицы является деревом и называется деревом исключения (elimination tree). Каждый узел дерева соответствует столбцу матрицы. Таким образом, мультифронтальный метод может быть представлен как обход дерева исключения от листьев к корню. При посещении очередного узла происходит построение фронтальной матрицы, в результате частичной факторизации которой алгоритм формирует соответствующий столбец фактора. Для построения фронтальной матрицы используются значения соответствующего столбца исходной матрицы, а также матриц обновления, ассоциированных с детьми рассматриваемого узла в дереве исключения. После выполнения частичной факторизации фронтальной матрицы формируется столбец фактора и матрица обновления, которая будет использована при построении фронтальной матрицы родителя. Высокоуровневое описание мультифронтального метода приведено ниже (алгоритм 1). Символом ® обозначена операция расширяющего сложения (extend add) [5].
Процедуры init_frontal_matrix, assembly_ frontal_matrix и form_update_matrix реализуются с помощью вызовов соответствующих функций BLAS. Более подробное описание метода можно найти в работах [5, 11]. Последовательная реализация мультифронального метода в решателе ННГУ описана в работе [12].
Недостатком базовой версии численной факторизации является неэффективное использование кэш-памяти. Для решения этой пробле-
мы на практике используются так называемые суперноды (supernode). Супернодом называется группа столбцов фактора, имеющих одинаковый шаблон ниже плотной треугольной подматрицы. Супернодальный мультифронтальный подход впервые был предложен Эшкрафтом, Гримсом, Льюисом, Пейтоном и Симоном [13] в 1987 г., а затем исследован Нг и Пейтоном [14]. Суперноды позволяют формировать фактор по несколько столбцов за итерацию с использованием функций BLAS третьего уровня, что повышает эффективность использования кэш-памяти. Именно эта редакция мультифрон-тального метода используется в качестве базовой для данной работы.
Алгоритм 1. Высокоуровневое описание муль-тифронтального метода
1 foreach node i of elimination tree
in topological order
2 init_frontal_matrix(^i)
3 foreach son j of i do
4 и ^ и eUj
5 end for
6 assembly _frontal_matrix(F, U)
7 factorize(F)
8 form_update_matrix(Ui)
9 Li^F(1*}
10 end for
СХЕМЫ РАСПАРАЛЛЕЛИВАНИЯ
В мультифронтальном методе могут быть использованы два способа организации параллелизма: применение параллельных функций BLAS и параллельное решение независимых друг от друга задач в соответствии со структурой дерева исключения. Рассмотрим перспективы применения этих методов.
Использование параллельного BLAS. Большая часть вычислений в мультифронталь-ном методе приходится на процедуры BLAS, такие как умножение плотных матриц и решение плотных систем линейных уравнений с треугольной матрицей. Поэтому использование существующих библиотек для высокопроизводительных вычислений, таких как, например, Intel MKL, является самым естественным способом распараллеливания численной фазы разложения Холецкого. К сожалению, эксперименты показывают, что применение этого подхода чаще всего приводит к разочаровывающим результатам. Этот факт объясняется тем, что большинство вспомогательных матриц, возникающих в ходе выполнения алгоритма, имеют маленькую размерность, и поэтому накладные
расходы, связанные с организацией параллелизма, не компенсируются последующей параллельной обработкой. Таким образом, необходимую производительность можно получить, только если в качестве основного метода распараллеливания использовать распараллеливание по дереву исключения.
Распараллеливание по дереву исключения. Распараллеливание мультифронтального метода может быть выполнено на основе дерева исключения Т, содержащего информацию о зависимостях по данным, возникающих в ходе вычислений. Пусть Т[к] - часть дерева исключения с корнем в вершине к. Показано [11], что два столбца 1 и ] фактора Ь могут быть вычислены параллельно тогда и только тогда, когда поддеревья Т[1] и Т[]] не пересекаются, т.е. не имеют общих узлов. Этот результат является основной для построения алгоритмов параллельной факторизации. При этом в качестве единицы работы (задачи) можно рассматривать построение фронтальной матрицы, соответствующей очередному узлу дерева. Основная проблема, возникающая при разработке параллельной редакции метода, заключается в наличии существенного дисбаланса между вычислительной трудоемкостью возникающих задач. Решение данных задач предполагает выполнение операций над матрицами, размеры которых могут отличаться на порядки от узла к узлу при сохранении общей тенденции укрупнения матриц при перемещении от листьев к корню. Для решения проблемы балансировки нагрузки могут использоваться как статические, так и динамические схемы.
Статические схемы распараллеливания. Существует множество методов, использующих статическую схему распараллеливания, однако большинство из них были разработаны для систем с распределенной памятью. В последнее время предпринимаются усилия [15] для переноса одного из наиболее эффективных методов статического распараллеливания - алгоритма Гейста-Нг [16] для работы в системах с общей памятью. Идея алгоритма заключается в нахождении некоторого слоя в дереве исключения, т.е. множества узлов, которые не обязательно находятся на одном уровне, но при этом не имеют общих потомков. Найденный слой должен обладать свойством сбалансированности, так чтобы количество операций, необходимых для обработки узлов поддеревьев с корнями в узлах найденного слоя, удовлетворяло заданному порогу и было примерно одинаковым. При реализации алгоритма свойство сбалансированности можно понимать как число операций
с плавающей точкой, либо как оценку времени, необходимого для обработки узла.
На рис. 1 приведен пример работы алгоритма. Найденный слой представляет срез дерева в узлах 6, 11, 14. Выделенные цветом поддеревья могут быть обработаны параллельно.
Рис. 1. Пример разбиения дерева исключения с использованием алгоритма Гейста-НГ. Разными цветами отмечены поддеревья, которые могут быть обработаны параллельно
двухуровневый параллельный
АЛГОРИТМ
Одним из основных недостатков статических схем является невозможность достаточно точно оценить объем работы, необходимый для обработки каждого узла дерева исключения. Поэтому в работе предлагается другой способ балансировки нагрузки, основанный на динамической схеме (рис. 2).
В рамках динамической схемы строится пул задач и на каждом шаге алгоритма поток достает задачу из очереди и приступает к ее выпол-
нению. Каждая задача соответствует вычислению соответствующего столбца фактора и состоит из четырех подзадач: вычисление матрицы узла, вычисление фронтальной матрицы, формирование из фронтальной матрицы столбца фактора, вычисление матрицы обновления. Пул задач организуется в виде приоритетной очереди.
Для формирования очереди задач используется алгоритм, который учитывает различные характеристики узлов дерева для достижения лучшей балансировки [3]. Алгоритм обходит все узлы дерева в соответствии с топологической перестановкой, сформированной раннее, и на каждой итерации цикла добавляет рассматриваемый узел в очередь. Приоритет узла в очереди складывается из основного и второстепенного, где первый отвечает за правильный обход столбцов в параллельном мультифронтальном методе, а второй - за улучшение балансировки.
Основной приоритет равен количеству детей узла в дереве исключения, а второстепенный вычисляется как оценка трудоемкости решения соответствующих подзадач. Трудоемкость решения задачи оценивается как количество операций с плавающей точкой.
Представленная схема позволяет эффективно использовать вычислительные ресурсы на нижних уровнях дерева исключения, но при приближении к корню возможности для параллелизма по задачам становятся ограниченными, при этом размеры обрабатываемых матриц значительно увеличиваются. В этот момент целесообразно изменить схему распараллеливания таким образом, чтобы вычислительные потоки использовали внутри вызовов многопоточные реализации функций BLAS (алгоритм 2).
Рис. 2. Пример работы мультифронтального метода на основе двухуровневой схемы распараллеливания
Алгоритм 2. Параллельный мультифронтальный метод на основе двухуровневой схемы
1 procedure process_node(node i of elimination_tree)
2 init_frontal_matrix(^)
3 foreach son j of i do
4 и ^ и eUj
5 end for
6 assembly_frontal_matrix(F,U
7 factorize(F)
8 form_update_matrix(Ui)
9 Li^F(1*)
10 end procedure
11
12 procedure two-level_parallel_multifrontal
13 omp_set_num_threads(MAX_SYSTEM_THREADS);
14 blas_set_num_threads(1);
15
16 #pragma omp parallel
17 while(task_queue.hasTaskET())
18 #pragma omp critical (queue)
19 i^task_queue.get_task_with_highest_priority();
20 process_node(i)
21 #pragma omp critical (queue)
22 task_queue.increase_task_primary_priority(parent(i))
23 end while
24
25 while(task_queue.hasTask())
26 i^task_queue.get_task_with_highest_priority();
27 process_node(i)
28 task_queue.increase_task_primary_priority(parent(i))
29 end while
30 end procedure
РЕЗУЛЬТАТЫ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ
Тестовая инфраструктура. Результаты экспериментов получены с использованием узла кластера, содержащего два восьмиядерных процессора Intel Sandy Bridge E5-2660 2.2 GHz, 64GB RAM, работающего под управлением ОС Linux CentOS 6.4. Использовался компилятор Intel C++ Compiler и библиотека Intel MKL BLAS из пакета Intel® Parallel Studio XE 2013 SP1. Для проведения экспериментов были выбраны матрицы из коллекции университета Флориды [17]. Характеристики тестовых матриц представлены в табл. Все они являются симмет-
ричными положительно определенными. В качестве перестановок, уменьшающих заполнение фактора, использовался METIS [18], но также могут быть использованы другие переупорядо-чиватели [19, 20].
Использование параллельных библиотек BLAS. Наиболее простым способом распараллеливания мультифронтального метода является использование высокопроизводительных параллельных библиотек BLAS. Для этого не требуется изменять исходный код, необходимо лишь собрать его с соответствующей реализацией BLAS. На диаграмме (рис. 3; слева) показано ускорение численной фазы разложения Холец-кого при запуске с различным числом потоков BLAS.
Таблица.
Характеристики тестовых матриц
Название Порядок Число ненулевых Число ненулевых
матрицы элементов элементов в факторе
Pwtk 217 918 5 926 171 49 025 872
Msdoor 415 863 10 328 399 51 882 257
parabolic fem 525 825 2 100 225 25 571 376
tmt_sym 726 713 2 903 835 28 657 615
boneS10 914 898 28 191 660 266 173 272
Emilia 923 923 136 20 964 171 1 633 654 176
audkiw 1 943 695 39 297 171 1 225 571 121
bone010 M 986 703 12 437 739 363 650 592
bone010 986 703 36 326 514 1 076 191 560
ecology2 999 999 2 997 995 35 606 934
thermal2 1 228 045 4 904 179 50 293 930
StocF-1465 1 465 137 11 235 263 1 039 392 123
Hook 1498 1 498 023 31 207 734 1 507 528 290
Flan 1565 1 564 794 59 485 419 1 451 334 747
G3_circuit 1 585 478 4 623 152 90 397 858
Каждый многоугольник на рисунке соответствует запуску с определенным количеством потоков. По осям отложено время работы на соответствующей матрице. Для матрицы Hook_1498 указано время работы в 8 потоков (при работе в 16 потоков запуск завершается с ошибкой из-за нехватки памяти).
Исходя из результатов запусков версии с параллельной библиотекой BLAS (рис. 3; слева) можно сделать следующие выводы. Все тестовые матрицы можно разделить на 2 группы в зависимости от степени их заполненности. Для первой группы, где матрицы больше и плотность их выше (матрицы Flan_1565, Emilia923, audikwl, bone010, StocF-1465,
Hook_1498), использование параллельной библиотеки BLAS позволяет получить ускорение до 6 раз при запуске в 16 потоков. Для второй группы, где матрицы либо маленькие, либо сильно разреженные (матрицы G3_circuit, pwtk, msdoor, parabolic_fem, tmt_sym, bone S10, bone010_M, ecology2, thermal2), ускорения при увеличении числа потоков BLAS практически не наблюдается. Также стоит отметить, что MKL BLAS имеет встроенный контроль размера матриц и не производит параллельную обработку, если размер матрицы слишком маленький. Таким образом, отсутствие ускорения на матрицах из второй группы можно считать хорошим результатом, т.к.при отсутствии подобного кон-
Рис. 3. Сравнение ускорения численной фазы разложения Холецкого: при запуске с разным числом потоков BLAS (слева); при запуске в 16 потоков с использованием параллельного BLAS
и динамической схемы (справа)
троля за размером матрицы можно было бы наблюдать замедление.
В работе [3] мы предложили использовать динамическую схему распараллеливания муль-тифронтального метода. Сравнивая ускорение с использованием различных схем распараллеливания (рис. 3; справа), можно видеть, что для 4 из 6 матриц первой группы (Emilia923, audikwl, bone010, Hook_1498) использование параллельного BLAS дает лучшее ускорение в среднем в 1,7 раза. Тем не менее, для остальных матриц предпочтительнее использовать динамическую схему. Этот факт говорит о том, что комбинация указанных методов распараллеливания (алгоритм 2) может дать потенциально лучшее ускорение по сравнению с каждой схемой в отдельности, что в дальнейшем подтверждается вычислительными экспериментами.
Выбор момента переключения между схемами распараллеливания. Важнейшим элементом алгоритма, влияющим на время работы численной фазы при использовании двухуровневой схемы распараллеливания, является критерий переключения между схемами (по узлам дерева; в рамках одного узла при помощи BLAS). Анализ параллельных запусков динамической схемы распараллеливания, использующей параллелизм на уровне логических задач в сочетании с последовательным BLAS, показал, что основным фактором, ограничивающим итоговое ускорение, является недостаток свободных задач, начиная с некоторого момента подъема по дереву исключения. В связи с этим предлагается использовать в качестве критерия сравнение числа необработанных задач с некоторым пороговым значением, выбираемым экспериментально. После срабатывания данного критерия происходит переход к последовательной обработке узлов дерева с использованием параллельного BLAS (алгоритм 2).
Для изучения вопроса о выборе порогового значения были выбраны 4 матрицы, представляющие 4 возможных класса: «небольшие разраженные» (parabolic_fem), «небольшие плотные» (pwtk), «большие разреженные» (G3_Circuit) и «большие плотные» (audikw_1).
Результаты приведены на диаграмме (рис. 4). По горизонтальной оси отложено пороговое значение, а по вертикальной - время работы двухуровневого алгоритма, отнесенное к времени работы с использованием динамической схемы (в зависимости от матрицы). Во всех запусках использовалось 16 вычислительных ядер. Значения, меньшие единицы, соответствуют преимуществу двухуровневой схемы над обычной.
Для небольших более плотных матриц и больших сильно разреженных время работы
при изменении порога не изменяется, кроме того, оно практически совпадает с временем работы при распараллеливании с использованием динамической схемы, отношение времен колеблется около единицы. Для больших более плотных матриц наблюдается значительное ускорение относительно динамической схемы, причем оно наблюдается уже при небольших значениях параметра и в дальнейшем практически не меняется.
Этот факт говорит о том, что большая часть работы сосредоточена в небольшой окрестности корня дерева, где параллелизм по задачам ограничен, но есть ресурс для использования параллельных библиотек BLAS. Ниже этой окрестности фронтальных матриц больше, они имеют средний размер, и использование параллелизма по задачам и параллелизма BLAS дает схожий результат. Для матриц из последней группы ситуация выглядит противоположным образом. Размер фронтальных матриц настолько мал, что использование параллельного BLAS сразу же приводит к замедлению. Однако таким замедлением можно пожертвовать, т.к. абсолютное время работы составляет менее 1 с.
На диаграмме (рис. 5) приведено абсолютное время работы мультифронтального метода с использованием различных схем распараллеливания. Во всех запусках использовались 16 ядер. Значение порога переключения схем для всех матриц было выбранным одинаковым и равнялось 200. Можно видеть, что для всех представленных матриц, время работы на которых превышает 5 с, новый двухуровневый метод показывает лучшие результаты. Данный эффект проявляется наиболее явно на матрице bone010, где удалось получить ускорение в 3 раза.
Отдельного внимания заслуживает вопрос об объеме используемой памяти. Так, фронтальные матрицы на верхних уровнях дерева исключения имеют больший размер. В связи с этим, при использовании параллелизма на задачах для обработки узлов верхних уровней требуется хранение вспомогательных структур данных для каждого из 16 потоков, что приводит к большим затратам памяти. Напротив, в двухуровневой схеме фронтальные матрицы узлов верхних уровней обрабатываются потоками совместно, что позволяет значительно сократить размер вспомогательных структур данных. В частности, применение описанных методов позволило получить правильное решение для матрицы Hook_1498 при использовании 16 потоков, в отличие от базовой динамической схемы.
Пример работы алгоритма. Проиллюстрируем работу алгоритма на примере (рис. 6). На рисунке изображены первые 11 уровней дерева
Рис. 4. Зависимость времени работы двухуровневого алгоритма от момента переключения между двумя параллельными схемами
Рис. 5. Сравнение времени работы мультифронтального метода при использовании динамической и двухуровневой схем
исключения в радиальной форме для матрицы StocF-1465. Корень дерева находится в центре, а концентрические круги соответствуют различным уровням дерева. Каждому узлу дерева поставлен в соответствие определенный цвет в зависимости от того, каким потоком он обрабатывался. Обратим внимание на то, что при приближении к корню дерева произошло переключение схемы работы на использование параллельного BLAS в сочетании с последовательной обработкой дерева, как и было заявлено в алгоритме 2. Представляет интерес распределение цветов, соответствующих разным потокам. Так, из рисунка видно, что в некоторых случаях в поддеревьях присутствуют цепочки узлов, отмеченных одним цветом. Такие примеры соответствуют достаточно эффективной схеме работы с памятью, допускающей повтор-
ное использование данных, ранее загруженных в кэш-память. Вместе с тем, наличие «разноцветных» поддеревьев показывает возможность дальнейшей оптимизации работы с памятью.
Сравнение с другими решателями. Был проведен ряд экспериментов на тестовых матрицах с использованием следующих известных и широко распространенных решателей:
• MKL PARDISO из Intel Math Kernel Library в составе Intel Parallel Studio XE 2013 SP1 [21];
• MUMPS (лучшее время работы из двух актуальных версий ver. 4.10.0, ver 5.0.0) [22].
Для всех пакетов использовались одинаковые перестановки, полученные с помощью METIS, а также функции BLAS и ScaLAPACK из библиотеки Intel MKL. Полученные результаты приведены на диаграмме (рис. 7).
Рис. 6. Пример работы двухуровневой схемы распараллеливания для матрицы 8йсР-1465
Из результатов можно сделать следующие выводы. Для запусков в 1 поток реализованный мультифронтальный метод показывает схожие результаты с МИМР8. Это объясняется тем фактом, что в обоих решателях в качестве метода численного разложения используется мультифронтальный метод. В то же время оба решателя проигрывают РАКОКО на 6 матрицах (ш8^ог, 1ш1_8уш, есо^у2, Шеташ12, 03_снсш1;, parabo1ic_fem), на 4 матрицах (аи&к^_1, Ьопе010, 8ЮсР-1465, Иап_1565) заметен выигрыш, в остальных случаях время работы сопоставимо. Наибольший выигрыш разработанного решателя по сравнению с РАКОКО получен на матрице аи&к—_1 и составляет 14 %, а наибольшее отставание, в 1,5 раза, на матрице есо^у2.
Для запусков в 4 потока на большинстве матриц наблюдается заметное отставание МиМР8 по сравнению с другими решателями. При этом разработанная двухуровневая динамическая схема распараллеливания показывает лучшие результаты по сравнению с РАКОКО на 6 матрицах (parabo1ic_fem, аи&к—_1, Ьопе010_М, Ьопе010, 8госБ-1465, Иап_1565), на остальных матрицах время работы обоих решателей сопоставимо.
При запуске в 8 потоков время работы РАКОКО и разработанного решателя практически для всех матриц сопоставимо, за исключением р-^к, ш8^ог, parabo1ic_feш и Ьопе810.
Однако, абсолютное время работы на этих матрицах небольшое, что приводит к недостаточной масштабируемости при использовании двухуровневой схемы. Решатель МИМР8 при запуске в 8 потоков значительно уступает двум другим решателям, особенно для запусков с большими матрицами (аи&к'_1, Иап_1565).
Наибольшее значение имеют результаты экспериментов при использовании 16 потоков. Можно видеть, что РАКО18О показывает преимущество во времени работы на 5 матрицах. Сравнивая двухуровневый алгоритм и МИМР8, нужно отметить преимущество первого на 6 матрицах (parabo1ic_feш, аи&к—_1,
Ьопе010_М, bone010, 8ЮсР-1465, Иап-1565). В свою очередь МИМР8 также выигрывает на 6 матрицах ф'-к, ш8^ог, 1ш1_8уш, eco1ogy2, 1Ьегша12, G3_circuit). Однако, если обратиться к диаграмме (см. рис. 5), можно видеть, что время работы на этих матрицах достаточно маленькое по сравнению с матрицами, на которых получен выигрыш.
Сравнивая время работы двухуровневого алгоритма и РАКОКО, можно отметить, что ситуация во многом схожая: матрицы, на которых отставание наиболее заметно ф-'к, ш8^ог, есо^у2, Шета^, G3_circuit), обрабатываются численно фазой достаточно быстро и по указанным ранее причинам дают худшее масштабирование. В остальных запусках время работы чис-
G3 circuit
Flan 1565
StocF-1465
therm а12
ecology2
1 Поток
pwtk
msdoor
G3 circuit
parabolicfem Flan_1565
tmtsym StocF-1465
boneSlO therm al2
': 1
audkiw 1
ecology2
4 Потока
pwtk
L6 ■
msdoor
parabolicfem
tmtsym
boneSlO
audkiw
boneOlO
boneOlO M
boneOlO
boneOlO M
8 Потоков
16 Потоков
pwtk
G3 circuit
msdoor
Flail 1565
StocF-1465
\ 1
therm al2
G3 circuit
parabolic_fem Flan_1565
tmt_sym StocF-1465
boneSlO therm al2
ecology2
audkiw 1
ecology2
msdoor
parabolic_fem
tmtsym
boneS10
audkiw 1
boneOlO
boneOlO M
boneOlO
boneOlO M
—PARDISO —MUMPS —Решатель, двухуровневая схема
Рис. 7. Сравнение численных фаз решателей. По осям отложено время работы. За единицу взято время работы MKL PARDISO
ленной фазы обоих решателей в большей степени сопоставимо.
ЗАКЛЮЧЕНИЕ
Основным результатом работы является комбинированная схема распараллеливания мультифронтального метода. Данная схема сочетает лучшие свойства двух подходов к организации параллелизма при обработке дерева исключения. Так, большое число «легковесных задач», соответствующих нижним уровням дерева, решается в рамках парадигмы распаралле-
ливания по задачам с динамической балансировкой нагрузки. При этом на верхних узлах дерева малое число «тяжеловесных задач» решается путем применения многопоточного BLAS. В работе сформулирован критерий переключения между схемами, показан выигрыш в производительности и памяти по сравнению с ранее подготовленной реализацией, выполнено сравнение с MKL PARDISO и MUMPS.
В дальнейшем планируется рассмотреть возможность усовершенствования разработанных программных реализаций за счет сокраще-
ния накладных расходов на организацию параллелизма. В частности, представляют интерес перспективы применения разных реализаций приоритетных очередей. Другим направлением дальнейших исследованием является разработка гетерогенных реализаций, использующих для расчетов не только традиционные процессоры, но и ускорители вычислений.
СПИСОК ЛИТЕРАТУРЫ
1. Li X. Direct Solvers for Sparse Matrices. [Электронный ресурс]. URL: http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ SparseDirectSurvey.pdf. [ X. Li. Direct Solvers for Sparse Matrices [Online]. Available: http://crd-legacy.lbl.gov/~xiaoye/ SuperLU/SparseDirectSurvey.pdf. ]
2. Козинов Е. А., Лебедев И. Г., Лебедев С. А., Мало-ва А. Ю., Мееров И. Б., Сысоев А. В., Филиппенко С. C. Новый решатель для алгебраических систем разреженных линейных уравнений с симметричной положительно определенной матрицей // Вестник Нижегородского университета им. Н. И. Лобачевского. 2012. № 5 (2). C. 376-384. [ E. A. Kozinov, I. G. Lebedev, S. A. Lebedev, et al., "Novel linear equations systems solver for sparse symmetric positive-definite matrix," (in Russian), in Vestnik UNN, no. 5 (2), pp. 376-384, 2012. ]
3. Lebedev S., Akhmedzhanov D., Kozinov E., Meyerov I, Pirova A., Sysoyev A. Dynamic Parallelization Strategies for Multifrontal Sparse Cholesky Factorization // LNCS, Parallel Computing Technologies. - Springer International Publishing, 2015. - pp. 68-79. [S. Lebedev, D. Akhmedzhanov, E. Kozinov, I. Meyerov, A. Pirova, A. Sysoyev. Dynamic Parallelization Strategies for Multifrontal Sparse Cholesky Factorization // LNCS, Parallel Computing Technologies. - Springer International Publishing, 2015. - pp. 68-79. ]
4. Davis Timothy A. Direct methods for sparse linear systems. Vol. 2. Siam, 2006. [T.A. Davis. "Direct methods for sparse linear systems", Vol. 2. Siam, 2006. ]
5. Liu J. W. The multifrontal method for sparse matrix solution: Theory and practice // SIAM review. 1992. Vol. 34, No. 1. P. 82-109. [J.W. Liu. "The multifrontal method for sparse matrix solution: Theory and practice", SIAM review. 1992. Vol. 34, No. 1. p. 82-109. ]
6. Kalinkin A., Arturov K. Asynchronous approach to memory management in sparse multifrontal methods on multiprocessors // Applied Mathematics. 2013. Vol. 4, No. 12A. P. 33-39. [ A. Kalinkin, K. Arturov. "Asynchronous approach to memory management in sparse multifrontal methods on multiprocessors", Applied Mathematics. 2013. Vol. 4, No. 12A. pp. 33-39. ]
7. George T., et al. Multifrontal factorization of sparse SPD matrices on GPUs // Parallel & Distributed Processing Symposium (IPDPS), 2011. P. 372-383. [ T. George, et al. "Multifrontal factorization of sparse SPD matrices on GPUs", Parallel & Distributed Processing Sympo-sium (IPDPS), 2011. p. 372-383. ]
8. Duff I.S., Reid J.K. The multifrontal solution of indefinite sparse symmetric linear // ACM Transactions on Mathematical Software (TOMS). 1983. Vol. 9, No. 3. P. 302-325. [ I.S. Duff, J.K. Reid. "The multifrontal solution of indefinite sparse symmetric linear", ACM Transactions on Mathematical Software (TOMS). 1983. Vol. 9, No. 3. pp. 302-325. ]
9. Liu J. W. The multifrontal method and paging in sparse Cholesky factorization // ACM Transactions on Mathematical Software (TOMS). 1989. Vol. 15, No. 4. P. 310-325. [ J. W. Liu. "The multifrontal method and paging in sparse Cholesky factorization", ACM Transactions on Mathematical Software (TOMS). 1989. Vol. 15, No. 4. pp. 310-325. ]
10. Amestoy P. R., et al. Vectorization of a multiprocessor multifrontal code // International Journal of High Performance Computing Applications. 1989. Vol 3, No. 3. P. 41-59. [P. R Amestoy, et al. "Vectorization of a multiprocessor multifrontal code", International Journal of High Perfor-mance Computing Applications. 1989. Vol 3, No. 3. pp. 41-59. ]
11. L'Excellent J. Y. Multifrontal Methods: Parallelism, Memory Usage and Numerical Aspects // Ph.D. thesis, Ecole normale superieure de lyon-ENS LYON. 2012. [ J. Y. L'Excellent "Multifrontal Methods: Parallelism, Memory Usage and Numerical Aspects", Ph.D. thesis, Ecole normale superieure de lyon-ENS LYON. 2012. ]
12. Лебедев С. А., Козинов Е. А. Разработка нового решателя разреженных систем линейных уравнений // Высокопроизводительные параллельные вычисления на кластерных системах: Материалы XIII Всероссийской конференции. / Изд-во ННГУ: Н. Новгород, 2013. С. 301-306. [ S. A. Lebedev, E.A. Kozinov. "Development of the new sparse SLAE solver", in Proc XIII International conference "High Performance computing on cluster systems". 2013. pp. 301-306. ]
13. Ashcraft C. C., et al. Progress in sparse matrix methods for large linear systems on vector supercomputers // International Journal of High Performance Computing Applications. 1987. Vol 1, No. 4. Pp. 10-30. [ C. C. Ashcraft, et al. "Progress in sparse matrix methods for large linear systems on vector supercomputers", International Journal of High Performance Computing Applications. 1987. Vol 1, No. 4. pp. 10-30. ]
14. Ng E., Peyton B.W. A supernodal Cholesky factorization algorithm for shared-memory multiprocessors // SIAM Journal on Scientific Computing. 1993. Vol 14, No. 4. P. 761769. [E. Ng, B. W. Peyton. "A supernodal Cholesky factorization algorithm for shared-memory multiprocessors", SIAM Journal on Scientific Computing. 1993. Vol 14, No. 4. pp. 761769. ]
15. L'Excellent J.Y., Sid-Lakhdar M.W. Introduction of shared-memory parallelism in a distributed-memory multifrontal solver // Research Report. 2013. RR-8227 <hal-00786055>. [ J.Y. L'Excellent, M.W. Sid-Lakhdar. "Introduction of shared-memory parallelism in a distributed-memory multifrontal solver", Research Report. 2013. RR-8227 <hal-00786055>. ]
16. Geist G., Ng E. Task scheduling for parallel sparse Cholesky factorization // International Journal of Parallel Programming. 1989. Vol. 18, No. 4. P. 291-314. [ G. Geist, E. Ng. "Task scheduling for parallel sparse Cholesky factorization", International Journal of Parallel Programming. 1989. Vol. 18, No. 4. pp. 291-314. ]
17. Davis T.A., Hu Y. The university of Florida sparse matrix collection // ACM Transactions on Mathematical Software (TOMS). 2011. Vol 38, No. 1. [ T.A. Davis, Y. Hu. "The university of Florida sparse matrix collection", ACM Transactions on Mathematical Software (TOMS). 2011. Vol 38, No. 1. ]
18. Karypis G., et al. A Fast and Highly Quality Multilevel Scheme for Partitioning Irregular Graphs // SIAM Journal on Scientific Computing. 1999. Vol. 20, No. 1. pp. 359-392. [G. Karypis, et al. "A Fast and Highly Quality Multilevel Scheme
for Partitioning Irregular Graphs", SIAM Journal on Scientific Computing. 1999. Vol. 20, No. 1. pp. 359-392. ]
19. Pellegrini F. Scotch and libScotch 6.0 User's Guide // Tech. rep., LaBRI. 2012. [F. Pellegrini. "Scotch and libScotch 6.0 User's Guide", Tech. rep., LaBRI. 2012. ]
20. Pirova A., Meyerov I. MORSy - a new tool for sparse matrix reordering // In Proceedings of an International Conference on Engineering and Applied Sciences Optimization. 2014. P. 1952-1963. [ A. Pirova, I. Meyerov. "MORSy - a new tool for sparse matrix reordering", In Proceedings of an International Conference on Engineering and Applied Sciences Optimization. 2014. pp. 1952-1963. ]
21. Intel Math Kernel Library Reference Manual [Электронный ресурс]. URL: https://software.intel.com/en-us/node/ 470282. [ Intel Math Kernel Library Reference Manual [Online]. Available: https://software.intel.com/en-us/node/ 470282 ]
22. MUltifrontal Massively Parallel Solver (MUMPS 5.0.0) User's guide [Электронный ресурс]. URL: http://mumps. enseeiht.fr/doc/userguide_5.0.0.pdf. [ MUltifrontal Massively Parallel Solver (MUMPS 5.0.0) User's guide [Online]. Available: http://mumps. enseeiht.fr/doc/userguide_5.0.0.pdf ]
ОБ АВТОРАХ
ЛЕБЕДЕВ Сергей Александрович, асп. каф. математического обеспечения и суперкомпьютерных технологий. М-р в обл. прикл. мат. и информатики (ННГУ, 2014). Иссл. в обл. разработки высокопроизводительных решателей разреженных СЛАУ.
МЕЕРОВ Иосиф Борисович, доц. каф. математического обеспечения и суперкомпьютерных технологий. М-р в обл. прикл. мат. и информатики (ННГУ, 1999). Канд. техн. наук по мат. моделированию, числ. методам и комплексам программ (ННГУ, 2005), доцент. Иссл. в обл. методов использования суперкомпьютеров для решения прикладных задач.
КОЗИНОВ Евгений Александрович, асс. каф. математического обеспечения и суперкомпьютерных технологий. М-р в обл. прикл. мат. и информатики (ННГУ, 2007). Иссл. в обл. разработки визуальных сред для параллельного программирования.
АХМЕДЖАНОВ Дмитирий Ринатович, студент каф. математического обеспечения и суперкомпьютерных технологий. Иссл. в обл. разработки высокопроизводительных решателей разреженных СЛАУ.
ПИРОВА Анна Юрьевна асп. каф. математического обеспечения и суперкомпьютерных технологий. М-р в обл. прикл. мат. и информатики (ННГУ, 2013). Иссл. в области методов переупорядочения разреженных матриц.
СЫСОЕВ Александр Владимирович, ст. преп. каф. математического обеспечения и суперкомпьютерных технологий. М-р в обл. прикл. мат. и информатики (ННГУ, 1999). Канд. техн. наук по мат. моделированию, числ. методам и комплексам программ (ННГУ, 2012). Иссл. в обл. параллельного программирования и глобальной оптимизации.
METADATA
Title: Two-level parallel strategy for multifrontal sparse Cholesky factorization.
Authors: S. A. Lebedev1, I. B. Meyerov2, E. A. Kozinov3, D. R. Akhmedzhanov4, A. Yu. Pirova5, A. V. Sysoyev6.
Affiliation:
1-6 Lobachevsky State University of Nizhni Novgorod (UNN), Russia.
Email: 1sergey.a.lebedev@gmail.com, 2meerov@vmk.unn.ru, 3evgeniy.kozinov@gmail.com, 4dmitry.akhmedzhanov@yandex.ru, 5anna.yu.malova@gmail.com, 6sysoyev@vmk.unn.ru.
Language: Russian.
Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 19, no. 3 (69), pp. 178-189, 2015. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print).
Abstract: In this paper we consider the problem of paralleliza-tion of Cholesky factorization numerical phase for sparse symmetric positive definite matrices. A new strategy for parallelization of the multifrontal method for shared-memory systems is suggested. This strategy combines two approaches to parallelism organization depending on the elimination tree level. At the bottom of the tree, parallel computing of nodes from a priority queue takes place. At the top levels of the tree, nodes are calculated sequentially, employing multithreaded BLAS procedures. Experimental results on the matrices from the University of Florida Sparse Matrix Collection are given. We show that our implementation is commensurable with MUMPS and MKL PARDISO solvers.
Key words: sparse algebra, Cholesky factorization, numerical phase, multifrontal method, high performance computing, dynamic parallelization, task-based parallelism.
About authors:
LEBEDEV, Sergey Alexandrovich, Post-graduate student of Software Department (UNN). Research in highperformance solvers of sparse SLAE. MEYEROV, Iosif Borisovich, PhD., associate professor of Software Department (UNN). Research in applications of supercomputers.
KOZINOV, Evgeny Alexandrovich, assistant teacher of Software Department (UNN). Research in visual environments for parallel programming.
AKHMEDZHANOV Dmitry Rinatovich, student of Software Department (UNN). Research in high-performance solvers of sparse SLAE.
PIROVA, Anna Yurievna, postgraduate student of Software Department (UNN). Research in ordering algorithms for the sparse matrices.
SYSOYEV, Alexander Vladimirovich, PhD., senior lecturer of Software Department (UNN). Research in parallel programming and global optimization.
