Системы линейных алгебраических уравнений (СЛАУ) с разреженной симметричной положительно определенной матрицей возникают при моделировании процессов во многих предметных областях. Огромная размерность таких систем (в современных приложениях -10 7 и выше) приводит к большим затратам памяти и процессорного времени, что без сомнения позволяет отнести их решение к области применения высокопроизводительных вычислений. Прямые методы, основанные на факторизации матрицы системы, активно применяются для решения больших разреженных СЛАУ. В настоящее время в этой области разработан целый ряд параллельных алгоритмов для систем с различной архитектурой и соответствующие программные пакеты [1], среди которых широко распространены MKL PARDISO, MUMPS, SuperLU, CHOLMOD и другие.

С 2011 года на факультете ВМК ННГУ разрабатывается прямой решатель разреженных СЛАУ с симметричной положительно определенной матрицей, основанный на методе Холецкого [2]. В рамках данного решателя изучаются вопросы оптимизации соответствующих алгоритмов под современные многоядерные архитектуры. В данной работе рассматривается задача распараллеливания наиболее затратной по времени и памяти численной фазы разложения Холецкого. Несмотря на наличие множества алгоритмов и их реализаций, вопрос о развитии существующих методов построения масштабируемых алгоритмов и программных средств, ориентированных на системы с общей памятью, не потерял свою актуальность в связи с постоянным развитием многоядерных архитектур. Ранее в работе [3] мы предложили способ распараллеливания мультифронтального метода, основанный на динамической схеме балансировки нагрузки. В данной статье предлагается модификация, состоящая в сочетании возможностей базовой схемы с использованием параллельного BLAS для решения «тяжеловесных» задач на верхних уровнях дерева исключения. Будет показано, что применение данной комбинированной схемы позволяет улучшить масштабируемость программной реализации.

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

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

Существует несколько методов выполнения численной фазы разложения Холецкого. Среди них можно выделить три наиболее широко используемых на практике: ориентированный влево (left-looking) [4], ориентированный вправо (right-looking) [4] и мультифронтальный (multifrontal) [5]. Основное отличие между ними заключается в способе формирования результирующей матрицы , а также в способе хранения и размещении в памяти промежуточных результатов. Перечисленные методы показывают схожую производительность на различных тестовых наборах матриц, но с точки зрения многих исследователей мультифронтальный метод является наиболее перспективным для распараллеливания [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 третьего уровня, что повышает эффективность использования кэш-памяти. Именно эта редакция мультифронтального метода используется в качестве базовой для данной работы.

end for

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

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

Распараллеливание мультифронтального метода может быть выполнено на основе дерева исключения T, содержащего информацию о зависимостях по данным, возникающих в ходе вычислений. Пусть T[k] -часть дерева исключения с корнем в вершине k. Показано [11], что два столбца i и j фактора L могут быть вычислены параллельно тогда и только тогда, когда поддеревья T[i] и T[j] не пересекаются, то есть не имеют общих узлов. Этот результат является основной для построения алгоритмов параллельной факторизации. При этом в качестве единицы работы (задачи) можно рассматривать построение фронтальной матрицы, соответствующей очередному узлу дерева. Основная проблема, возникающая при разработке параллельной редакции метода, заключается в наличии существенного дисбаланса между вычислительной трудоемкостью возникающих задач. Решение данных задач предполагает выполнение операций над матрицами, размеры которых могут отличаться на порядки от узла к узлу при сохранении общей тенденции укрупнения матриц при перемещении от листьев к корню. Для решения про-блемы балансировки нагрузки могут использоваться как статические, так и динамические схемы.

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

На рисунке 1 приведен пример работы алгоритма. Найденный слой представляет срез дерева в узлах 6, 11, 14. Выделенные цветом поддеревья могут быть обработаны параллельно. Для формирования очереди задач используется алгоритм, который учитывает различные характеристики узлов дерева для достижения лучшей балансировки [3]. Алгоритм обходит все узлы дерева в соответствии с топологической перестановкой, сформированной раннее, и на каждой итерации цикла добавляет рассматриваемый узел в очередь. Приоритет узла в очереди складывается из основного и второстепенного, где первый отвечает за правильный обход столбцов в параллельном мультифронтальном методе, а второй -за улучшение балансировки.

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

Представленная схема позволяет эффективно использовать вычислительные ресурсы на нижних уровнях дерева исключения, но при приближении к корню возможности для параллелизма по задачам становятся ограниченными, при этом размеры обрабатываемых матриц значительно увеличиваются. В этот момент целесообразно изменить схему распараллеливания таким образом, чтобы вычислительные потоки использовались внутри вызовов многопоточные реализации функций BLAS (алгоритм 2). [17]. Характеристики тестовых матриц представлены ниже (таблица 1). Все они являются симметричными положительно определенными. В качестве перестановок, уменьшающих заполнение фактора, использовался METIS [18], но также могут быть использованы другие переупорядочиватели [19,20]. В работе [3] мы предложили использовать динамическую схему распараллеливания мультифронтального метода. Сравнивая ускорение с использованием различных схем распараллеливания (рисунок 3; справа) можно видеть, что для 4 из 6 матриц первой группы (Emilia923, audikw1, bone010, Hook_1498) использование параллельного BLAS дает лучшее ускорение в среднем в 1.7 раза. Тем не менее, для остальных матриц предпочтительнее использовать динамическую схему. Этот факт говорит о том, что комбинация указанных методов распараллеливания (алгоритм 2) может дать потенциально лучшее ускорение по сравнению с каждой схемой в отдельности, что в дальнейшем подтверждается вычислительными экспериментами.

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

Для изучения вопроса о выборе порогового значения были выбраны 4 матрицы, представляющие 4 возможных класса: «небольшие разраженные» (parabolic_fem), «небольшие плотные» (pwtk), «большие разреженные» (G3_Circuit) и «большие плотные» (audikw_1).

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

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

Был проведен ряд экспериментов на тестовых матрицах с использованием следующих известных и широко распространенных решателей:

 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. Полученные результаты приведены на диаграмме (рисунок 6). Рассматривая результаты экспериментов для запусков решателей в 16 потоков можно видеть, что PARDISO сохраняет преимущество во времени работы и показывает лучшие результаты на 5 матрицах. Сравнивая двухуровневый алгоритм и MUMPS, нужно отметить преимущество первого на 6 матрицах (parabolic_fem, audikw_1, bone010_M, bone010, StocF-1465, Flan-1565). В свою очередь MUMPS также выигрывает на 6 матрицах (pwtk, msdoor, tmt_sym, ecology2, thermal2, G3_circuit). Однако, если обратиться к диаграмме (рисунок 5), можно видеть, что время работы на этих матрицах достаточно маленькое по сравнению с матрицами, на которых получен выигрыш. Сравнивая время работы двухуровневого алгоритма и PARDISO можно отметить, что ситуация во многом схожая: матрицы, на которых отставание наиболее заметно (pwtk, msdoor, ecology2, thermal2, G3_circuit), обрабатываются численно фазой достаточно быстро и по причинам, описанным в предыдущем разделе, дают худшее масштабирование. В остальных запусках время работы численной фазы обоих решателей в большей степени сопоставимо.

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