CC BY
102
УДК 519.612
0 НЕКОТОРЫХ ВАРИАНТАХ МЕТОДА
___ _________ ________________ _________«_»
ДЕКОМПОЗИЦИИ ОБЛАСТЕЙ1 В.П. Ильин, Д.В. Перевозкин
Рассматриваются алгоритмы масштабируемого распараллеливания решения сверхбольших разреженных сеточных СЛАУ, представленных в универсальных сжатых форматах, в том смысле, что их реализация осуществляется без программных ограничений на порядки алгебраических систем и на количество используемых вычислительных узлов, процессоров и/или ядер.
Данная задача сводится к распределенному варианту алгебраической 3Б-декомпозиции областей, в котором отсутствует чрезмерная расчетно-информационная нагрузка корневого процессора, т.е. все организуемые MPI-процессы, каждый из которых соответствует своей подобласти, являются практически равноправными. Вычислительный процесс состоит из двух основных этапов, первый из которых заключается в непосредственной автоматической декомпозиции, на основе анализа матричного портрета и формировании крупноблочного представления СЛАУ. Второй этап — это реализация крыловского итерационного алгоритма FGMRES (гибкого обобщенного метода минимальных невязок), использующего точное или приближенное обращение диагональных матричных блоков (многопоточное решение подсистем в подобластях с использованием средств OpenMP) с помощью прямого или итерационного метода соответственно.
Описываемые методы реализованы в составе библиотеки алгебраических решателей Krylov. В работе приводятся некоторые оценки используемых ресурсов и особенности параллельных вычислительных технологий. Эффективность разработанных алгоритмов иллюстрируется результатами численных экспериментов по решению характерных алгебраических задач на различных конфигурациях многопроцессорной вычислительной системы.
Ключевые слова: декомпозиция областей, матричный граф, распараллеливание алгоритмов, сеточные уравнения, разреженные СЛАУ, программные и вычислительные технологии
Введение
Методы декомпозиции областей (МДО), как главное орудие масштабируемого распараллеливания численного решения больших задач математического моделирования на многопроцессорных вычислительных системах (МВС), является предметом изучения огромного количества статей и монографий, представленных, например, на специальном интернет-сайте [1].
Обширную область изучаемых здесь вопросов можно разбить на две основные части. Первая из них относится к построению новых алгоритмов для различных классов задач, их теоретическому обоснованию, оценкам скорости сходимости, оптимизации и сравнительному анализу. Относительно свежие математические результаты в данных направлениях можно найти на сайте 22-й Международной конференции по декомпозиции областей [2], состоявшейся в сентябре 2013 г.
Вторая часть не менее актуальных проблем — это эффективное отображение МДО на архитектуру современных МВС, которая, несмотря на стремительное обновление суперком-пьютерного списка ТОР-500, в ближайшую пятилетку, до рождения первого экзафлопсни-
1 Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2014»
ка, прогнозируется традиционной: кластер из гетерогенных узлов, содержащих по несколько центральных многоядерных процессоров (СРи) с общей памятью, а также по несколько графических процессоров общего назначения (ОРОРИ), сверхбыстрых, но с достаточно сложным использованием. Управление иерархической памятью с распределенным и общим доступом ограничивается средствами МР1, ОрепМР и СИБЛ.
Цель данной работы заключается в том, чтобы хотя для достаточно ограниченного набора алгоритмов декомпозиции сконцентрироваться на вариантах технологических решений, обеспечивающих реализацию в определенном смысле главной задачи — создания высокопроизводительного математического и программного обеспечения, не содержащего формальных ограничений на число степеней свободы и на количество используемых вычислительных ядер и/или процессоров.
Рассматриваемая нами задача с функциональной точки зрения достаточно простая — это численное решение систем линейных алгебраических уравнений (СЛАУ) с разреженными матрицами больших порядков (до 108—1011), возникающими при различных аппроксимациях многомерных междисциплинарных начально-краевых задач со сложной геометрической конфигурацией границ на неструктурированных сетках.
Особенность проблемы заключается в объеме данных, необходимых для хранения ненулевых элементов разреженных матриц, число которых в каждой из строк не зависит от размерности СЛАУ и составляет обычно несколько десятков (их количество растет с увеличением порядка точности аппроксимации). Например, система уравнений Навье —Стокса, описывающая гидро-газодинамические течения, включает 4 неизвестные скалярные функции — три компоненты вектора скорости и давление. Поэтому на трехмерной сетке с числом узлов 10003 = 109 общий порядок СЛАУ составляет примерно 4 ■ 109, а количество ненулевых матричных элементов — около 1011. Поскольку возникающие при этом вычислительные задачи являются очень плохо обусловленными, то расчеты необходимо проводить, как минимум, со стандартной двойной точностью (64 двоичных разряда, или 8 байтов, на одно число). Таким образом, размещение соответствующей СЛАУ на одном процессоре, объем оперативной памяти которого составляет обычно несколько гигабайт, становится проблематичным. В этом случае не проходит зачастую применяемая стратегия распараллеливания алгебраических задач, когда сначала глобальная матрица формируется на одном (головном) процессоре, а потом ее блоки рассылаются по остальным используемым процессам. Единственная альтернатива здесь — распределенное хранение матрицы, когда ее элементы изначально рассчитываются в разных процессорах.
Данная статья построена следующим образом. В разделе 1 мы представляем основные структурные особенности рассматриваемых СЛАУ и применяемых МДО. В разделе 2 описываются некоторые предлагаемые параллельные алгоритмы и их программные реализации. В разделе 3 приводятся результаты характерных численных экспериментов по решению представительной серии методических задач.
1. Некоторые постановки, формализмы и алгоритмы декомпозиции
Рассматривается следующая задача: пусть матрица СЛАУ
Аи = /, А = [аг,3 }еП*’*, (1)
разбита на блочные строки
А = {Ар е Пи?'и, р = 1,Р}, N1 + ... + = N, (2)
с приблизительно равным числом строк N ^ 1 в каждой. Соответствующим образом также
разбиваются векторы искомого решения и правой части
и = {ир еП^}, / = {/р еПм*}, (3)
так что система уравнений (1) может быть записана в форме совокупности Р подсистем
р
Ар,рир + Ху Ар,дия = /р, р = 1, •••, Р, (4)
9 = 1
9 = Р
где АР;д е '%,МР’Мя — прямоугольные матричные блоки, получаемые при разбиении каждой блочной строки (и всей матрицы А) на блочные столбцы:
Ар = {Ар,д; Ч = 1,...,Р}, А = {Ар,д;р, Ч = 1,...,Р}. (5)
Будем считать, что СЛАУ (1) представляет собой систему сеточных уравнений, аппрок-
симирующих некоторую многомерную краевую задачу для дифференциального уравнения, так что каждая компонента векторов и, / соответствует узлу сетки, общее число которых
р
в расчетной сеточной области = У 0р равно N. Отметим, что в методах конечных
р=1
элементов (МКЭ, [4]) такие системы называются узловыми, но более сложные сеточные уравнения мы пока не рассматриваем. При этом блочное разбиение матриц и векторов соответствует разбиению (декомпозиции) на Р сеточных непересекающихся подобластей 0р, в каждой из которых находится Жр узлов, N1 + ... + Nр = N.
Описанная декомпозиция не использует узлов-разделителей, т.е. условные границы подобластей не проходят через узлы сетки. Формальности ради можно считать, что внешность расчетной области есть неограниченная подобласть 0о с числом узлов N0 = 0. Сеточное уравнение для г-го узла сетки может быть записано в виде
аг,ги; + I] а;,,и,- = /г, г е 0й, (6)
З=г
где через Ш; обозначается сеточный шаблон г-го узла, т.е. совокупность номеров всех узлов (как правило, являющихся геометрическими соседями), участвующих в г-м уравнении.
Если сеточные узлы и соответствующие переменные пронумерованы следующим образом: сначала идут подряд все узлы 1-й подобласти О1 (неважно в каком внутреннем порядке), затем — узлы второй подобласти и т.д., — то вышеприведенное блочное разбиение матриц и векторов, приводящее к подсистемам (4), соответствует декомпозиции расчетной области без пересечения подобластей. При этом номера узлов из Ш;, соседних к г-му узлу, указываются в СБИ-формате [3] для г-й строки матрицы (фактически там указываются номера столбцов как для диагонального, так и внедиагональных элементов). Отметим, что взаимнооднозначного соответствия алгебраической и геометрической (сеточной) интерпретации структуры СЛАУ может и не существовать. Типичный пример — одному узлу
сетки может соответствовать т > 1 переменных в алгебраической системе. Соответству-щая ситуация имеет место в междисциплинарных задачах, описываемых системами из т дифференциальных уравнений, у которых решение есть вектор-функция размерности т.
При заданном блочном разбиении СЛАУ, или декомпозиции сеточной расчетной области без пересечения подобластей, для решения системы (1) можно построить хорошо распараллеливаемый аддитивный метод Шварца, представимый итерационным процессом
ар,риП = /р- £ Ар,чиП-1 = #р-1
9=1 9=Р
При вычисленных правых частях ^рЩ-1 нахождение всех иЩ требует для каждой подобласти решения подсистем с порядками ^, которые на каждой п-й итерации могут выполняться независимо и одновременно на разных процессорах. В данном случае перевычисление д^-1 фактически означает использование новых значений иЩ-1 для соседних подобластей в качестве граничного условия 1-го рода (Дирихле) для околограничных внутренних узлов из 0р.
Рассмотренное двойственное представление структуры СЛАУ (1) может быть дополнено путем введения графового описания этой же задачи. Каждому сеточному узлу с номером г (или г-й строке матрицы А) можно сопоставить вершину V; некоторого графа С, а матричному элементу аг,з — ребро графа С = (V, Е) ив соответствии с [5] определить множества вершин и ребер:
V = {V;; г = 1,N}, Е = {(^,и,-Ж, = 0}. (8)
Отметим следующий существенный технологический аспект, возникающих при реализации численного решения подсистемы (7). Традиционно универсальные программные решатели СЛАУ, ориентированные на матричный СБИ-формат, предполагают сплошную нумерацию векторных компонент и соответствующих матричных строк, начиная с единицы или нуля. Однако, если для исходной СЛАУ информация дана в глобальной нумерации переменных, то в подобласти 0р требуется своя локальная нумерация. Таким образом, здесь возникает техническая задача переформирования С£Ер-формата для диагонального блока Ар,р с переходом из глобальной нумерации в локальную. А после вычисления подвекто-ров ир в подобластях, естественно, потребуется сборка глобального вектора и, т.е. перевод векторных компонент из локальной нумерации в глобальную.
Известно, что скорость сходимости итераций метода Шварца (7) возрастет, если декомпозицию расчетной области сделать с пересечением подобластей.
В силу этого мы дадим определение расширенной сеточной подобласти 0р Э 0р с пересечениями, величину которых мы будем определять в терминах количества околограничных сеточных слоев, или фронтов.
Обозначим через Г0 е 0р множество внутренних околограничных узлов из 0р, т.е. таких узлов Р; е 0р, у которых один из соседей не лежит в 0р(Р, е 0р, ] е Шг, ] = г)). В Гр определим подмножество узлов Гр,д, у которых соседние по шаблону узлы принадлежат смежной подобласти 0д, ч е 0°, где 0 — совокупность номеров подобластей, соседних с 0р. Таким образом,
Г0 = II Г0
х р ^ х р,д>
причем подмножества Г°,^ могут пересекаться, т.е. содержать околограничные узлы, у которых есть соседи из различных смежных подобластей.
Обозначим далее через Гр множество узлов, соседних с узлами из Гр, но не принадлежащих 0р, через Гр — множество узлов, соседних с узлами из Гр, но не принадлежащих объединению Гр и 0р, через Г3 — множество узлов, соседних с Гр, но не принадлежащих
(7)
Гри Гри 0Р, и т.д. Соответственно эти множества назовем первым внешним слоем (фронтом) узлов, вторым слоем, третьим и т.д. Получаемое объединение узлов
= 0р и гр... и Г^ (10)
будем называть расширенной р-й сеточной подобластью, а целую величину А определяем
как величину расширения, или пересечения (в терминах количества сеточных слоев). Случай А = 0 фактически означает декомпозицию области 0^ на подобласти без пересечений (00 = 0р).
Множество ГА е 0^ представляет совокупность внутренних околограничных узлов расширенной подобласти 0^, а Г^+1 есть множество внешних околограничных узлов. Условная граница 0^ проходит между Г^ и Г^+1.
Аналогично Гр, множество Г^ можно разбить на подмножества околограничных узлов
ГА = ГА и ГА и ГА (11)
Р Р,9тр \±±1
у которых соседние (по шаблону) находятся соответственно в подобластях 091, 0?2,..., 09тр (здесь тр означает количество подобластей, пересекающихся с 0А, а 91, 92,..., 9тр — номера этих подобластей). Отметим, что ГА может содержать «кратные» околограничные узлы, принадлежащие одновременно разным подмножествам Гр,д, т.е. у которых есть соседи из различных смежных подобластей.
Если обозначить через 0р^д пересечение какой-то подобласти 0д с расширенной подобластью 0р, то последнюю можно представить в форме объединения
0р = 0р и 0Р,,1 ••• и 0А9т1, (12)
причем каждое пересечение составляется из частей внешних околограничных фронтов:
°Ад = ГР,9 и ГР,9... и ГМ, 9 = . (13)
Таким образом, сборка матрицы СЛАУ для расширенной подобласти может быть интерпретирована как идентификация, объединение и перенумерация множества узлов, участвующих в (12), (13).
На формальном алгебраическом языке расширенные подсистемы уравнений и итерационный метод Шварца записываются аналогично (7):
Ар,рир = /Р — ^ Ар,дир 1 = $р Р. (14)
9=1
9 =Р
Здесь размерность векторов и /р равна числу узлов N в расширенной подобласти Ор. При этом надо иметь в виду, что в пересечениях один узел принадлежит двум или более подобластям.
Остановимся теперь подробнее на формировании правых частей др1-1 или д^-1 в (7),
(13). Для этого рассмотрим сеточное уравнение вида (4) для расположенного в Ор і-го околограничного узла, у которого некоторые из соседей находятся в других подобластях Од ,<? = р:
аг,гйг + ^ аг,іЙ = /і — X} Й• (15)
Отметим, что в (15) г-й узел располагается и в Пр, и в какой-то из соседних подобластей. Отсюда строим итерационный процесс, несколько модифицировав данное уравнение:
(а»,г + О Е агЛ < + Е аг,.?и™ = /г + Е (О<-1 - и™-1). (16)
у ^'/Пр 7 ^'епр ^-(/пр
Здесь О £ [0,1] — некоторый итерационный параметр, который при О = 0 соответствует условию Дирихле в околограничном г-м узле, при О = 1 — условию Неймана, а при
0 < О < 1 — условию 3-го рода, или Робена. Значения и™-1 , и™-1 в правой части (16) берутся с предыдущей итерации в узлах из соседних к Пр подобластей. Очевидно, что если итерационный процесс (16) сходится, то при любом О его предел будет решением исходной системы (15).
Если ввести диагональную ^р = ^гад{ Е аг,^}, то итерационный процесс (16) можно
3 = ^р
записать в форме, обобщающей соотношение (14):
Вр(иП - ип-1) = г^-1 = /р - (Аип-1 )р, (17)
где предобуславливающая матрица Вр имеет вид Вр = Ар,р + О^р.
В приведенных выше рассуждениях мы рассмотрели расширение подобласти Пр в наружную сторону. Однако такие же процедуры осуществляются и для соседних подобластей, что приводит к построению пересечений внутри самой Пр. Их реализация может осуществляться по сеточных слоям (фронтам), как и при расширении Пр.
Остановимся коротко на вопросах декомпозиции при решении нестационарных задач, когда использование неявных аппроксимаций по времени приводит к необходимости решения СЛАУ на каждом шаге. В этом случае приходится решать последовательно алгебраические системы с разными правыми частями и одинаковыми матрицами, если их элементы не меняются со временем.
2. Особенности реализации распределенных алгоритмов
В высокопроизводительных версиях декомпозиции областей существуют два ключевых момента: формирование разреженных алгебраических систем для соответствующих подобластей и реализация параллельных по подобластям крыловских итерационных процессов якобиевого типа.
2.1. Метод построения квазисбалансированных сеточных подобластей
В начале данного раздела остановимся на одном факторе, имеющем большое значение для скорости сходимости итерационных процессов вида (7) по подобластям. Речь идет
о топологии и некоторых других характеристиках разбиения расчетной сеточной области П^- на сеточные подобласти П^-, р = 1,...,Р. Здесь следует отметить, что вообще говоря, нужно различать исходную непрерывную расчетную область П, в которой решаются дифференциальные и/или интегральные уравнения, и ее дискретный аналог П^. Более того, декомпозицию можно делать в первую очередь: сначала разбить область на подобласти
( П = У Пр ), а затем в каждой из них строить подсетку П^, которая геометрически привя-^ р=1 7
зана к соответствующей Пр. Однако мы эти аспекты сознательно опускаем, ограничиваясь только декомпозицией сеточной области.
Совокупность подобластей образует собой макросеть, или макрограф, у которого каждая подобласть представляет собой макровершину, а связи с геометрически соседними подобластями — макроребра. Как известно [6], скорость сходимости итераций по подобластям будет тем выше, чем меньше диаметр макрографа (максимальное из минимально возможных расстояний между макроузлами). При разбиении трехмерной области на достаточно большое число Р подобластей диаметр макрографа есть 0(Р) при одномерной декомпозиции, 0(Р1/2) — при двумерной и 0(Р1/3) — при трехмерной. Поэтому естественно, что оптимальной является трехмерная сбалансированная декомпозиция, при которой число узлов в каждой подобласти является примерно одинаковым, чтобы при одновременном решении алгебраических подсистем в подобластях соответствующие процессоры имели минимальные простои.
В некотором смысле альтернативным подходом к декомпозиции сеточной области является следующий. В соответствии с вышеприведенными рассуждениями, оптимальное разбиение кубической области с кубической же сеткой — это декомпозиция на кубические сеточные подобласти. Особенность этой конфигурации заключается в том, что при заданном числе узлов она имеет относительно малый диаметр соответствующего сеточного подграфа, т.е. сравнительно близкий к идеальному варианту — сфере. Именно это наблюдение мы и положим в основу рассматриваемого «локального» принципа декомпозиции: строить сеточные подобласти с примерно равным количеством узлов и имеющие относительно малые диаметры своих сеточных подграфов.
Итак, описываемый алгоритм декомпозиции преследует цель создания сеточных подобластей Пр, р = 1, ...,Р, по возможности меньшего диаметра, где диаметр понимается как диаметр подграфа Ор(рр,Ер), порожденного [7] множествами содержащихся в подобласти вершин рр и ребер Ер. Вообще говоря, поиск оптимального в терминах какого-либо критерия разбиения графа О является задачей как минимум полиномиальной сложности, хотя для решения больших СЛАУ нежелательна даже квадратичная зависимость от количества элементов графа Nу, NЕ.
С этой точки зрения интересен следующий ход: с помощью какого-либо экономичного метода выполнить на первом этапе разбиение графа на связные непересекающиеся подмножества сравнительно небольшой мощности и малого диаметра, впоследствии рассматривая их как макровершины и декларируя наличие ребра между двумя соседними макровершинами в том и только том случае, если есть хотя бы одно ребро, соединяющее некоторые элементы этих двух подмножеств. В дальнейшем предлагается интерпретировать полученный таким образом граф как огрубленную (агрегированную) копию исходного, применяя более сложные (или, наоборот, простые) алгоритмы для формирования последующих и конечного разбиений.
Не останавливаясь пока на методах агрегации каждого этапа, отметим, что предлагаемый подход является многошаговым и фактически образует последовательность вложенных сеток, или макросеток, из подобластей. Если обозначить через в номер шага в рассматриваемой иерархической структуре, то для формализации данного процесса можно ввести следующие обозначения: О(в)(У(в), Е(в)), в = 0,1,..., М, — граф макросетки в-го уровня (в = 0 соответствует начальному сеточному графу, в котором вершина представляет простой узел сетки, а М есть количество стадий агрегации, которое заранее неизвестно), ^ и — число макровершин и макроребер макросетки в-го уровня.
В целях построения начального разбиения предлагается использовать поиск в ширину [5]. Данный метод обладает привлекательным свойством помечать в первую очередь вершины с наименьшим расстоянием от исходной, что соответствует цели создания подмножеств малого диаметра. Поиск ведется среди вершин, еще не относящихся к какому-либо подмножеству, до тех пор, пока количество помеченных вершин не достигнет определенного предела (желаемого размера макровершины N^ax, где Nhax означает число узлов сетки, которое допускается включать в макровершину первого уровня). Все помеченные на j-м шаге алгоритма вершины объявляются принадлежащими j-му подмножеству, после чего процедура повторяется до исчерпания множества вершин.
Формализуем алгоритм поиска в ширину. Через Adj(v) обозначим множество вершин, смежных вершине v. Пусть также Q — структура данных типа очередь [5]. Обозначим через
Q ^ S добавление множества S к очереди (в произвольном порядке), v ^ Q — изъятие
элемента из очереди. Для выполнения поиска в ширину дополнительно введем C и W — целочисленные массивы (списки) цветов и весов вершин (смысл последних будет раскрыт позднее). Через C(v) и W(v) обозначим текущие цвет и вес вершины v. Перед стартом алгоритма цвет всех вершин инициализируется нулем, а их вес полагаем равным единице. Также введем nmax — желаемый размер макровершины. Тогда алгоритм заключается в начальном выборе произвольной вершины v из множества V ив дальнейших преобразованиях графа в соответствии с псевдокодом, представленным на рис.
i = 1
while {u £ V | C(u) = 0} = 0 pick any v from {u £ V | C(u) = 0}
Q := v
n=0
while (n < nmax and Q = 0) v ^ Q C(v) := i
Q ^ Adj(v) П {u £ V | C(u) = 0} n = n + W (v) end while
i = i + 1 end while
После работы алгоритма объявим V = {u | C(u) = i}. Сформируем граф с макровершинами G1 = (V;, E;), где V; = {Vi,... Vn'}, N — количество получившихся макровершин, а множество ребер E определяется так:
E' = {(Vi, Vj) | 3u £ Vi, v £ Vj | (u, v) £ E} . (18)
Указанный выше алгоритм можно использовать и для формирования собственно подобластей нужного размера, положив W(Vi) = |Vi|, а nmax равным желаемому размеру подобласти. На данный момент именно такой алгоритм используется в библиотеке Krylov, хотя, вообще говоря, представляется целесообразным использование некоторого аналога очереди с приоритетом, где, например, приоритет тем выше, чем больше ребер связывают рассматриваемую макровершину с текущим составом формируемой подобласти, и тем ниже,
чем меньше вес (мощность) этой макровершины (стоит заметить, что в таком случае разумно вычислять приоритет макровершин в момент извлечения элемента из очереди, а не в момент добавления). Такой алгоритм будет отдавать предпочтение макровершинам с размером меньше типичного, выполняя их поглощение, и тем макровершинам, включение которых способствует уменьшению «площади» границы подобласти (т.е. количества ребер, связывающих эту подобласть с остальными).
Предполагается, что созданные таким образом подобласти будут обеспечивать хорошую скорость сходимости метода Шварца.
2.2. Особенности программной реализации
Рассмотрим более подробно некоторые вопросы, касающиеся непосредственно организации вычислительного процесса, причем остановимся для краткости на алгоритмах без пересечения подобластей.
В первую очередь заметим, что выполнение одного шага аддитивного метода Шварца (14) состоит из P независимых операций решения СЛАУ, для реализации каждой из которых формируется MPI-процесс, а выполнение ее осуществляется одним из процессоров CPU с указываемым количеством доступных вычислительных ядер, или потоков Nth, которыми выполняется внутреннее распараллеливание над общей памятью средствами OpenMP.
Во вторых, реально итерации по подобластям проводятся не по алгоритму Шварца
(14), а с помощью предобусловленного гибкого обобщенного метода минимальных невязок FGMRES [8], в котором блочно-диагональная матрица с блоками Ap,p на главной диагонали выступает в качестве предобуславливателя. Решение подсистемы для каждой p-й подобласти осуществляется прямым алгоритмом или итерационным, и в последнем случае предобуславливатель фактически является переменным.
Учитывая все вышесказанное, логичной представляется следующая схема реализации вычислений. Будем хранить векторы так же, как и матрицу, по Np компонент в каждом процессе, соответствующем своей подобласти. В таком случае векторно-векторные операции и обращение предобуславливающего оператора можно выполнять без обмена данными между MPI-процессами (исключение составляет лишь скалярное произведение, где объем пересылаемых данных невелик). Для вычисления действия оператора A на вектор будем пересылать лишь те данные, которые требуются для вычисления соответствующих частей вектора, т.е. только между соседними подобластями.
Итерации по подобластям продолжаются до выполнения условий
IHI = II/ - Au”1l < ^ 1 или n < nmax, (19)
где ее и nmax — задаваемые критерии останова. Если решение вспомогательных СЛАУ в подобластях осуществляется итерационным методом, то для него задаются внутренние критерии ei, nmax.
3. Примеры численных экспериментов
Главная задача данного раздела заключается в том, чтобы продемонстрировать эффективность итерационного метода по подобластям, которая в значительной степени зависит от качества декомпозиции. Расчеты проводились на кластере ИВМиМГ СО РАН [9] для
количества подобластей P = 1, 2, 4, 8,16, 32, причем случай P = 1 соответствует решению задачи без применения декомпозиции.
В табл. 1 мы приводим результаты по решению модельной задачи Дирихле для диффузионно-конвективного уравнения в декартовой системе координат ж, y, z с постоянными конвективными коэффициентами p = q = r = 16, см. [10, 11], в кубической области Q = [0,1]3, аппроксимируемой на кубической сетке с числом узлов N = N с помощью экспоненциальной семиточечной схемы. Граничные условия соответствуют точному решению и(ж, y, z) = 1, а начальное приближение выбрано и0(ж, y, z) = 0. Для критерия останова (19) внешнего итерационного процесса во всех случаях было выбрано ee = 10-7.
Отметим, что в приведенных расчетах использовался двухуровневый вариант декомпозиции (s = 0,1 в разделе 2), причем выбирались значения №2Х = vN, = N/P.
В каждой клетке табл. 1 представлены сверху вниз: количество внешних итераций, астрономическое и процессорное время (ta и tcpu, последнее — суммарное по всем MPI-процессам и по всем вычислительным потокам). В качестве решателя в подобластях использовался PARDISO из библиотеки Intel MKL [3], при заданном числе потоков Nth = 12. Из-за высоких требований данного решателя к объему оперативной памяти оказалось невозможным получить на имеющемся оборудовании результаты, соответствующие незаполненным клеткам таблицы.
Таблица 1
Решение модельной задачи с использованием PARDISO в подобластях
N \ P 1 2 4 B 1б 32
643 1 32,2 1б4 29 14,2 93,5 32 б,39 55,1 B3 39 ,B б, 324 54 3,42 232 7B 1,58 200
І283 1 BB5 10340 40 221 49B7 44 Вб,б 17B5 53 30,1 7B2 75 20,4 1232 10B 12,б 1б15
2563 — — б0 2222 100009 72 332 13450 102 212 1B400 142 13B 1B210
В табл. 2 приводятся результаты для решения той же задачи, однако для решения СЛАУ в подобластях использовался предобусловленный метод бисопряженных градиентов со стабилизацией (БЮ081аЪ), причем для предобуславливания использовалась модификация Айзенштата [4, 8] алгоритма неполной факторизации, без использования какого-либо распараллеливания. Критерии останова полагались следующими: = 10-1, п^ах = 10.
В табл. 3 представлена зависимость между критериями останова ег, птах, фактически отвечающих за точность решения СЛАУ в подобластях, и величинами, отражающими скорость сходимости внешнего итерационного процесса, а также общую вычислительную сложность. В каждой клетке таблицы перечислены: количество внешних итераций пе, общее количество внутренних итераций П (суммарное по всем подобластям и по всем внешним итерациям), астрономическое и процессорное время в секундах. Для сравнения также
Таблица 2
Решение модельной задачи с использованием предобусловленного итерационного метода БЮ081аЪ в подобластях, ег = 0,1, п^ах = 10
N \ Р 1 2 4 8 16 32
6 со 9 4.05 4.05 29 6,35 12,7 34 3,44 13,2 47 1,54 12,3 65 1,59 20,9 98 2,26 77,3
1283 18 64.9 64.9 58 115 230 56 57,4 229 67 25,2 201 92 15,9 215 132 18,3 510
2563 21 606 606 64 1045 2079 74 647 2588 92 443 3553 143 271 3916 92 231 6627
приводятся величины, получаемые при решении СЛАУ в подобластях прямым методом РЛИБШО.
Таблица 3
Влияние точности решения СЛАУ в подобластях на скорость сходимости и ресурсоемкость
Задача \ критерии ^ = 0,1 Птах 10 ^ = 0, 01 Птах = 20 ^ = 0,001 Птах = 40 РАИБШО
67 57 54 53
8 3297 5280 7040 424
N = 1283 25,2 39,9 52,2 30,1
201 318 415 782
143 111 108 102
Р = 16 15929 18843 25489 1632
N = 2563 271 319 427 212
3916 5112 6824 18400
Отметим, что во всех приведенных экспериментах итоговая погрешность численного решения была примерно одинаковой и равнялась 10-6, что достаточно хорошо соответствует критерию ее = 10-7.
По приведенным результатам можно сделать следующие выводы.
• С увеличением количества подобластей Р и соотвествующих им МР1-процессов, до определенного числа наблюдается существенное ускорение решения СЛАУ как для прямого, так и итерационного алгоритмов в подобластях. Однако, в силу роста числа внешних итераций пе с увеличением Р, а также ввиду роста количества обменов при этом, наблюдается минимум £а(Р), который на грубой сетке с N = 643 наступает при Р = 8, а с увеличением N проявляется при гораздо большем числе Р. Уменьшения числа пе можно достичь за счет применения декомпозиции с пересечением подобластей и использования других приемов ускорения итераций по подобластям, однако
вопрос о возможности монотонного роста ускорения вычислений с ростом P остается открытым.
• Сравнение времен решения СЛАУ при использовании для решения СЛАУ в подобластях прямого или итерационного алгоритмов показывает, что в различных ситуациях выигрывает то один, то другой подход. Как и следовало ожидать из теоретических оценок вычислительной сложности, прямой алгоритм проигрывает итерационному, если число узлов в подобласти достаточно велико. Однако с ростом значений P ресурсо-емкость прямого метода убывает быстрее, чем итерационного. Справедливости ради следует заметить, что времена счета соответствуют реализации прямого алгоритма с распараллеливанием по 12 потокам, а итерационного — без распараллеливания.
• Что касается решения СЛАУ с применением итерационного алгоритма в подобластях, то здесь существенные резервы ускорения вычислений заложены в оптимизации критериев окончания внутренних итераций. Как видно из табл. 3, уменьшение и увеличение n^ax приводит к излишним итерациям в подобластях и увеличению времени счета. Приведенные результаты носят предварительный характер и оставляют открытым вопрос об оптимальной стратегии выбора критериев останова внутренних итераций.
Заключение
Предложенный и экспериментально исследованный в работе алгоритм автоматизированного построения сбалансированной декомпозиции, как видно из расчетов, позволяет достаточно эффективно реализовывать параллельное решение СЛАУ на системах с распределенной памятью. Сравнение прямых и итерационных методов решения вспомогательных СЛАУ в подобластях свидетельствует о преимуществе итерационных процессов в подпространствах Крылова, особенно при решении сверхбольших алгебраических систем. При этом решение подсистем на каждой внешней итерации достаточно проводить с невысокой точностью.
Работа поддержана грантом РФФИ № 14-07-00128, а также грантами Президиума РАН № 15.9-4 и ОМН РАН № 1.3.3-4.
Литература
1. Domain Decomposition Methods. URL: http://ddm.org (дата обращения: 14.03.2012)
2. 22nd International Conference on Domain Decomposition Methods (DD22). URL: http: //dd22.ics.usi.ch/ (дата обращения: 31.11.2013)
3. Intel (R) Math Kernel Library from Intel. URL: http://software.intel.com/en-us/ articles/intel-mkl/ (дата обращения : 08.04.2014)
4. Ильин, В.П. Методы и технологии конечных элементов / В.П. Ильин — Новосибирск: Изд-во ИВМиМГ СО РАН, 2007. — 371 с.
5. Писсанецки, С. Технология разреженных матриц / С. Писсанецки — М.: Мир, 1988. — 305 с.
6. Bramble, J.H. Convergence estimates for product iterative methods with applications to domain decomposition / J. Bramble, J. Pasciak, J. Wang, J. Xu // Mathematics of
Computation. — 1991. — V. 57, № 195. — P. 1-21.
7. Берж К. Теория графов и ее применения / К. Берж — М.: Изд-во иностр. лит., 1962. — 250 c.
8. Saad, Y. Iterative Methods for Sparse Linear Systems, Second Edition / Y. Saad — SIAM, 2003. — 528 p.
9. Кластер НКС-30Т. URL: http://www2.sscc.ru/HKC-30T/HKC-30T.htm (дата обращения:
12.02.2013).
10. Andreeva, M.Yu. Two solvers for nonsymmetric SLAE / M.Yu. Andreeva, V.P. Il’in, E.A. Itskovich // Bulletin NCC, «Numerical Analysis» series. — 2003. — Iss. 12. — P. 1-16.
11. Бутюгин, Д.С. Методы параллельного решения СЛАУ на системах с распределенной памятью в библиотеке Krylov / Д.С. Бутюгин, В.П. Ильин, Д.В. Перевозкин // Вестник ЮУрГУ. Серия «Вычислительная математика и информатика». — 2012. — Т. 47, № 306. — С. 5-19.
Валерий Павлович Ильин, д.ф.-м.н., профессор, главный научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), ilin@sscc.ru.
Данил Валерьевич Перевозкин, младший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), foxillys@gmail.com.
Поступила в редакцию 10 апреля 2014 г.
Bulletin of the South Ural State University Series “Computational Mathematics and Software Engineering”
2014, vol. 3, no. 2, pp. 5-19
ON SOME VARIANTS OF DOMAIN DECOMPOSITION METHODS
V.P. Il’in, Institute of Computational Mathematics and Mathematical Geophysics SB
RAS (Novosibirsk, Russian Federation),
D.V. Perevozkin, Institute of Computational Mathematics and Mathematical
Geophysics SB RAS (Novosibirsk, Russian Federation)
The paper considers the algorithms for solving large sparse SLAEs arising from grid approximations of boundary value problems. The SLAEs and algorithms are not limited in a sense of number of unknowns, computational nodes, processors and/or cores. This problem is reduced to a distributed variant of algebraic 3D-domain decomposition, in which no excessive load of the root process is present, i.e. all MPI-processes, each of which corresponds to its own subdomain, are almost equal. The computational process consists of two main stages. The first stage is the automatic decomposition, based on the analysis of the matrix portrait and the formation of large-block representation of the original SLAE. The second stage implements a Krylov subspace iterative process with FGMRes (flexible generalized minimal residual method) using either exact or approximate inverse of diagonal blocks as a preconditioner. The methods described are implemented as a part of Krylov, a library of algebraic solvers. The paper presents some features of current parallel implementation and estimates of resource usage. Efficiency of the developed algorithms is illustrated by solving several typical model problems with different parameters and in different configurations of multiprocessor computer systems.
Keywords: domain decomposition methods, matrix graphs, parallel algorithms, grid equations, sparse linear algebraic systems
References
1. Domain Decomposition Methods. URL: http://ddm.org (accessed: 14.03.2012)
2. 22nd International Conference on Domain Decomposition Methods (DD22). URL: http: //dd22.ics.usi.ch/ (accessed: 31.11.2013)
3. Intel (R) Math Kernel Library from Intel. URL: http://software.intel.com/en-us/ articles/intel-mkl/ (accessed: 08.04.2014)
4. Il’in V.P. Metody i tekhnologii konechnykh elementov [Methods and technologies of finite elements]. Novosibirsk, ICM&MG SBRAS Publishing, 2007.
5. Pissanetski S. Tekhnologiya razrezhennykh matrits [Sparse matrix technology]. Moscow, Mir Publishing, 1988.
6. Bramble, J.H. Convergence estimates for product iterative methods with applications to domain decomposition / J. Bramble, J. Pasciak, J. Wang, J. Xu // Mathematics of Computation. — 1991. — V. 57, № 195. — P. 1-21.
7. Berzh K. Teoriya grafov i eye primeneniya [Graph theory and its applications]. Мoscow: Foreign Literature Publishing, 1962.
8. Saad, Y. Iterative Methods for Sparse Linear Systems, Second Edition / Y. Saad — SIAM, 2003. — 528 p.
9. NKS-30T cluster. URL: http://www2.sscc.ru/HKC-30T/HKC-30T.htm (accessed:
12.02.2013).
10. Andreeva, M.Yu. Two solvers for nonsymmetric SLAE / M.Yu. Andreeva, V.P. Il’in,
E.A. Itskovich // Bulletin NCC, «Numerical Analysis» series. — 2003. — Iss. 12. — P. 1-16.
11. Butyugin D.S., Ilin V.P., Perevozkin D.V. Metody parallelnogo resheniya SLAU na sistemakh s raspredelennoy pamyatju v biblioteke Krylov [Parallel SLAE solution methods for distributed memory systems in Krylov library] Vestnik Yuzho-Uralskogo gosudarstvennogo universiteta. Seriya "Vychislitelnaya matematika i informatika"[Bulletin of South Ural State University. Series: Computational Mathematics and Informatics]. 2012. No. 47(306). P. 5-19.
Received 10 April 2014-
