CC BY
178
УДК 519.63
МЕТОДЫ ПАРАЛЛЕЛЬНОГО РЕШЕНИЯ СЛАУ НА СИСТЕМАХ С РАСПРЕДЕЛЕННОЙ ПАМЯТЬЮ В БИБЛИОТЕКЕ KRYLOV1
Д.С. Бутюгин, В.П. Ильин, Д.В. Перевозкин
Рассматривается подход к созданию итерационного black-box («черного ящика») параллельного решателя, использованный в библиотеке Krylov для систем линейных алгебраических уравнений (СЛАУ) с разреженными матрицами высокого порядка, возникающими при сеточных аппроксимациях многомерных краевых задач и представленными в сжатом строчном формате CSR. Предлагается вариант алгебраической одномерной декомпозиции СЛАУ. Алгоритм основан на обходе в ширину графа матрицы системы и позволяет привести ее к блочно-трехдиагональному виду. За основу алгебраического решателя системы взят аддитивный метод Шварца, который естественным образом ложится на архитектуру вычислительных систем с распределенной памятью. Полученные алгебраические системы в подпространстве следов, образованных переменными на внутренних границах подобластей, решаются с помощью обобщенного метода минимальных невязок. Вспомогательные системы в подобластях решаются с помощью прямого алгоритма PARDISO из библиотеки Intel MKL, использующего распараллеливание над общей памятью средствами OpenMP. Реализованные алгоритмы апробированы на численном решении ряда задач вычислительной математики, таких как задачи гидродинамики, диффузионно-конвективные уравнения, задачи электромагнетизма и др. Приведенные результаты численных экспериментов демонстрируют эффективность предлагаемых решений для многопроцессорных вычислительных систем с распределенной памятью.
Ключевые слова: итерационные алгоритмы, методы декомпозиции областей, распараллеливание, алгебраические системы, разреженные матрицы, численные эксперименты, аддитивный метод Шварца.
Введение
Системы линейных алгебраических уравнений с разреженной матрицей возникают при численном решении большого числа задач вычислительной математики. Вычислительная сложность таких задач часто требует кластерных расчетов. В настоящее время существует ряд подходов к декомпозиции задач и, соответственно, матриц систем, а также итерационных алгоритмов, позволяющих свести решение исходной задачи к серии решений задач в подобластях. К методам декомпозиции задачи можно отнести геометрические методы декомпозиции расчетной области и последующую аппроксимацию задачи в каждой из подобластей, а также алгебраические методы декомпозиции матрицы системы, например, реализованные в библиотеке METIS [1].
В рамках данной работы предлагается вариант алгебраической одномерной декомпозиции СЛАУ. Алгоритм основан на обходе в ширину графа матрицы системы и позволяет привести ее к блочно-трехдиагональному виду. За основу алгебраического решателя системы взят аддитивный метод Шварца, который естественным образом ложится на архитектуру вычислительных систем с распределенной памятью. К данному методу применяется крыловское ускорение [2], заключающееся в замене стационарного процесса xn+ = Txn + g решением системы (I — T)x = g, которое
1 Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии 2012»
осуществляется при помощи итерационных методов в подпространствах Крылова, таких как GMRES. При этом итерации осуществляются в подпространстве следов, образованных переменными на границах подобластей. Решение систем в подобластях осуществляется параллельным прямым решателем разреженных систем PARDISO из библиотеки Intel® MKL [3]. Использование такого подхода позволяет реализовать black-box решатель, не требующий дополнительной информации для решения СЛАУ и имеющий небольшое число конфигурационных параметров. Реализация алгоритма на языке C+—+ с использованием средств MPI и высокоуровневой библиотеки линейной алгебры Eigen [4] включена в библиотеку Krylov [5].
Тестирование алгоритма проводилось на большом числе методических и практических задач вычислительной физики. Было произведено исследование быстродействия решателя и числа итераций в зависимости от различных параметрах решателя. К варьируемым параметрам относились параметры алгебраической декомпозиции, такие как перекрытие подобластей и их количество, а к параметрам времени исполнения — число потоков в PARDISO при решении задач в подобластях, а также число MPI процессов на узел. Кроме того, была продемонстрирована точность получаемых решений для представленных задач.
1. Алгоритмы решения СЛАУ
1.1. Метод алгебраической декомпозиции СЛАУ
Аппроксимация различных краевых задач соответствующих дифференциальных уравнений различными методами, такими как методы конечных разностей, конечных объемов и конечных элементов, приводит к системе линейных алгебраических уравнений вида
Ax = b, x,b Е Rn, A E Rn,n, det |A| = 0. (1)
Решение СЛАУ на кластере требует предварительной ее декомпозиции и распределения блоков матрицы, а также векторов неизвестных и правой части между узлами кластера. В работе предлагается использовать следующую модификацию метода одно-направленной алгебраической декомпозиции задачи [2], основными достоинствами которой являются простота, высокая скорость работы алгоритма и возможность получить распределенную блочно-трехдиагональную матрицу на выходе.
Сначала ищется псевдо-периферийная вершина v графа G, ассоциированного с матрицей системы (1).
Алгоритм 1. Поиск псевдо-периферийной вершины v графа G. Запускается обход в ширину [6] с произвольной вершины графа матрицы и находятся наиболее удаленные от нее вершины и расстояние до них. Далее снова запускается обход в ширину, начиная уже с произвольной вершины среди множества наиболее удаленных вершин, полученных на предыдущем шаге. Обходы в ширину повторяются до тех пор, пока увеличивается максимальное расстояние от стартовой вершины до остальных. Как только расстояние перестанет увеличиваться, алгоритм останавливается и в качестве псевдо-периферийной вершины v возвращает вершину, с которой начинался последний обход в ширину.
Время работы алгоритма 1 равно O(E ■ K), где E — количество ребер в графе матрицы. Для большинства конечных схем аппроксимации уравнений число ребер
пропорционально количеству вершин сетки в расчетной области. Этот факт следует из того, что соответствующие конечные базисы финитны и, как правило, используются сетки, у которых имеются ограничения на минимальные величины плоских и телесных внутренних углов. Число К в оценке равно количеству итераций, совершенных алгоритмом 1. Для произвольных графов в худшем случае это число может быть довольно велико, однако для графов практических задач алгоритму обычно требуется всего лишь несколько итераций [2]. За счет этого на практике алгоритм поиска псевдо-периферийных вершин оказывается весьма эффективным. Более того, во многих случаях оказывается достаточно выбрать произвольную, например, первую вершину графа, и вообще не запускать алгоритм поиска псевдо-периферийных вершин.
После того, как выбрана некоторая вершина V, начиная с нее запускается обход графа О в ширину. Все вершины разделяются на множества — «фронты» — согласно расстоянию от них до вершины V. То есть, во фронт (к = 0.. .в,, где в — псевдо-диаметр графа О) попадают все вершины, равноудаленные от вершины V на расстояние к.
Далее рассмотрим способ объединения фронтов в алгебраические подобласти Пр. Мы будем объединять только последовательные фронты в подобласти, так что
гр
^р = и ^к >
к—
где 1р и гр — границы подобласти Пр. Пусть требуется построить Р подобластей (например, по числу используемых узлов кластера) для системы (1) с N неизвестными. Желательно разбить неизвестные в подобласти приблизительно по N/P неизвестных. Однако поскольку мы объединяем в подобласти не отдельные вершины, а фронты, то точное разбиение неизвестных по N/P в каждой подобласти может не получиться. Тем не менее, мы можем при помощи бинарного поиска минимизировать максимальный размер подобласти либо максимизировать минимальный размер.
Рассмотрим, например, максимизацию минимального размера, которая может использоваться для того, чтобы, по возможности, более равномерно загрузить все узлы кластера. Легко видеть, что функция зависимости максимально возможного числа подобластей при фиксированном минимально допустимом размере монотонна. Действительно, при фиксированном минимальном размере разделение вершин графа разбиение на подобласти можно произвести при помощи следующего жадного [6] алгоритма.
Алгоритм 2. Алгоритм построения максимального числа подобластей. К первой подобласти добавляем последовательно фронты, начиная с первого, пока размер подобласти не станет больше либо равным минимально допустимому размеру. В этом момент можно прекратить увеличение первой подобласти и приступить к построению второй подобласти. Этот процесс продолжается до тех пор, пока все фронты не окажутся в какой-либо подобласти. При этом может оказаться так, что в конце окажется несколько фронтов, которых недостаточно для формирования отдельной подобласти. Такие фронты можно отнести к последней построенной подобласти.
Ясно, что приведенный выше способ строит максимальное число подобластей. Действительно, если бы можно было построить больше подобластей при заданном
минимальном размере подобласти, то тогда, начиная с некоторого г, подобласть П в «жадном» разбиении и подобласть П' отличались бы. Рассмотрим минимальное такое г. Если в подобласти П' больше фронтов, чем в П^, то тогда оптимальное разбиение можно перестроить таким образом, чтобы в подобласти г было столько же фронтов, что и в оптимальном. Значит, в подобласти П' должно быть меньше фронтов, чем в П^. Однако П строилось таким образом, что в него поместили минимально возможное число фронтов при фиксированном начальном. Получаем противоречие, и значит, построенное алгоритмом разбиение является максимальным.
Из приведенных выше рассуждений становится ясно, что с увеличением минимального размера число подобластей не может увеличиться. Значит, функция зависимости их числа от минимального размера действительно монотонна, и можно использовать двоичный поиск максимального минимального размера подобласти, чтобы получить Р подобластей. Кроме того, описанный выше жадный алгоритм позволяет собственно построить эти подобласти уже после определения минимального размера.
Представленный алгоритм генерирует подобласти без пересечений. Однако получить разбиение с пересечениями довольно легко — достаточно просто расширить каждую из подобластей на несколько фронтов слева и справа. Также можно при начале генерации следующей подобласти отступить от границы предыдущей внутрь необходимое число фронтов.
Заметим, что ребра в графе матрицы могут соединять вершины либо одного, либо двух соседних фронтов. Тогда получается, что граничными для подобласти р будут вершины, принадлежащие двум фронтам ^р-1 и ^Гр+1. Можно считать, что эти фронты принадлежат соседним подобластям Пр-1 и Пр+1 соответственно. Получается, что если последовательно занумеровать переменные в каждой из подобластей, а в каждой из подобластей — последовательно в каждом из фронтов, мы получим следующую блочно-трехдиагональную СЛАУ:
Ах =
В и Ь1 В2
Ьр-1
ир-1
Вр
х1 /1
х2 /2 = /
хр /р
Будем полагать, что диагональные блоки Вр, соответствующие СЛАУ в подобластях Пр, не вырождены. При этом если разбиение исходной СЛАУ выполнялось с пересечениями, то в полученной СЛАУ (2) порядка N/ окажется по несколько переменных, соответствующих одной и той же вершине графа, но лежащим в разных подобластях; в этом случае N1 > N.
1.2. Крыловское ускорение аддитивного метода Шварца
После декомпозиции задачи на подобласти и пересылки блоков матрицы в соответствующие узлы кластера мы имеем распределенную блочную СЛАУ (2). Для решения таких СЛАУ широко используются методы Шварца. Одна итерация аддитивного метода Шварца в матричном виде может быть записана как (см. [2])
Хі+1 = Xі + ^ ЯГВ-1Яр(/ - Ах*),
р
где Яр : ^ — оператор сужения на подобласть Пр. Здесь N1 — размерность
СЛАУ (2), а N — количество неизвестных в подобласти Пр. С использованием блочного вида СЛАУ итерацию метода Шварца можно переписать следующим образом:
жр+1 = жр + Бр 1(/р — ЬрЖр-1 — ирхр+1),
где для р =1 и р = Р полагаем, что Ьржр-1 = 0 и иржр+1 = 0 соответственно. В общем случае члены Ьржр-1 и /ржр+1 являются аналогами условий Дирихле для краевых задач.
Заметим, что на самом деле итерировать достаточно только переменные, образующие внешние границы подобластей, поскольку, зная граничные значения, легко определить значения внутри подобластей, решив соответствующую подзадачу. Для этого в произведениях Ьржр-1 и [/ржр+1 достаточно использовать только части векторов из подобластей Пр-1 и Пр+1, соответствующие фронтам 1р — 1 и гр + 1. Эти переменные и будут образовывать внешнюю границу области Пр.
Обозначим за пь вектор граничных переменных для всех подобластей, и введем оператор сужения на граничные переменные Я : 'Я.м ^ 'Я.Мь, где N — размерность пространства следов. Тогда итерация метода Шварца в подпространстве следов выглядит следующим образом:
пЬ+1 = пЬ + Яяр- с-1Яр(/ — пЬ).
р
Эту итерацию можно представить в виде пЬ+1 = Б(пь) = Тпь + д, где
т = I — Я£ярТДр-1ЯрЛяЯт, д = Я£яртДр_1Яр/. рр Сходимость стационарных итераций имеется при спектральном радиусе р(Т) < 1. В то время как для эллиптических уравнений это условие выполнено, для произвольных уравнений это не так. Однако можно заметить, что решение пь должно удовлетворять п = Тп* + д, что равносильно системе
[I — Т] пь = д. (3)
Данную систему можно решать итерационными методами в подпространствах Крылова, например обобщенным методом минимальных невязок (GMR.ES) [2]. Для этого достаточно, чтобы матрица [I — Т] = Я^ Яр^Д—1ЯрАЯт была невырожденной, что является существенно более слабым условием, по сравнению с условием р(Т) < 1. Заметим также, что матрица [I — Т] не является симметричной, поэтому нельзя использовать экономичные методы в подпространствах Крылова, такие как метод сопряженных градиентов, ориентированные на решение симметричных СЛАУ и использующие короткие рекурсии для построения базиса подпространств Крылова.
Метод GMRES рассчитан на решение произвольных невырожденных СЛАУ вида Аж = Ь, при этом он минимизирует невязку г, = Ь — Аж, в подпространстве Крылова
Кт(А, Го) = 8рап(го, Аго, А2го,..., Ат-1го}.
Псевдокод алгоритма следующий:
• Го — b — Axo, в = ||ro||, qo — ro/||ro||, С — (1,0,
• j — 0,1,... while ||rj|| > ев
- qj+1 — Aqj
- For k G [0,..., j]
* Hk,j — (qj+1,qk)
* qj+1 — qj+1 — Hk,j qk
- Hj+1,j — ||qj+1||, qj+1 — qj+1/Hj+1,j
- For к G [0,..., j — 1]
0)
T
* k — ck s k k
Hk+1,j sk ck Hk+1,j
cj — | j |/V |j | 2 + |Hj+1,j 2
kj — cjHj+1,j/Hj,j
0 1 cj sj 0
_ 0+1 —— —sj cj . 0+1
Hj,j — cj Hj,j
sj Hj+1,j, Hj+1,j — 0
- ||г^'+1|1 ^ в|С^'+1|
• у ^ вн-1е
• X ^ Хо + [до ... д,-1]у'.
Из псевдокода видно, что для алгоритма требуется только операция умножения матрицы системы на вектор. Запишем результат умножения I — Т на вектор V через результат одной итерации метода аддитивного Шварца Б(V):
[I — Т] V = V — ^ = V — Б^)+ д, д = Б(0).
Таким образом, для решения СЛАУ (3) методом GMR.ES достаточно реализовать обычный метод аддитивного Шварца. Для решения задач в подобластях хр = 1ЬР
можно использовать какой-нибудь прямой решатель для разреженных систем. При этом Ьи-разложение матриц Др достаточно вычислить всего один раз.
Оценим эффективность данного алгоритма для простого случая — решения трехмерного уравнения Пуассона в параллелепипеде со структурированной сеткой. Пусть имеется сетка с п вершинами по каждой из координатных осей. Тогда матрица системы будет иметь порядок N = п3. При решении такой системы прямыми методами требуется О^5/3) памяти и О^7/3) операций [7], поскольку ширина ленты матрицы к в этом случае оказывается пропорциональной п2.
Предположим теперь, что декомпозиция задачи на Р подобластей выполнена равномерно с небольшими перехлестами, то есть порядок СЛАУ в каждой подобласти есть О^/Р). Размерность пространства следов будем считать равной О(Р • п2) = О(Р • ^/3). Пусть для решения СЛАУ (3) требуется К итераций. Считая, что в подобластях сохраняются оценки для прямых решателей, получим, что общие затраты памяти М будут равны
М = О(Р • ^/Р)5/3 + К • Р • N2/3) = О^5/3/Р2/3 + К • Р • N2/3), или, в пересчете на один узел (при условии использования распределенного GMRES),
Мр = О(^/Р)5/3 + К • N2/3).
Время работы решателя Т будет складываться из максимального времени факторизации £fact, времени работы итерационного метода Ьцег и времени на коммуникации
tcomm (T = tfact + titer + tcomm). Для НИХ МОЖНО ПОЛуЧИТЬ Следующие Оценки:
tfact = O((N/P)7/3),
titer = O(K ■ (N/P)5/3 + K2 ■ N2/3),
tcomm = O(K ■ P ■ N2/3).
Данные оценки справедливы для распределенного GMRES. При использовании GMRES на одном узле оценка для titer увеличиваются до
titer = O(K ■ (N/P)5/3 + K2 ■ P ■ N2/3).
Тогда имеем
T = O((N/P)7/3 + K ■ (N/P)5/3 + K2 ■ P ■ N2/3).
Поскольку, как правило, K растет с ростом P, получаем, что при фиксированном N для T имеется оптимум при некотором значении P. Этот оптимум достигается, с одной стороны, за счет уменьшения времени факторизации блоков, и, с другой стороны, за счет увеличения временных затрат на ортогонализацию векторов и межпроцессные коммуникации.
2. Технологии реализация решателя
При выборе языка программирования, сторонних библиотек и технологий для реализации решателя учитывалась специфика используемых в решателе алгоритмов. Необходим был высокоуровневый язык и набор библиотек, хорошо подходящий для реализации различных алгоритмов на графах, манипуляций с векторами и разреженными матрицами. Дополнительное требование заключалось в том, чтобы компиляторы языка генерировали высоко-эффективный код. В итоге выбор был сделан в пользу языка C+—+ и использования стандартной библиотеки STL, из которой использовались различные контейнеры, а также библиотеки линейной алгебры Eigen [4]. Последняя имеет в своем составе классы разреженных матриц и динамических векторов, предоставляет удобный API для манипуляции ими, а также имеет встроенную векторизацию алгоритмов с использованием инструкций SSE.
Для решения задач в подобластях требуется эффективный решатель для СЛАУ с разреженными матрицами невысокого порядка, причем необходимо решать СЛАУ с различными правыми частями. Кроме того, обусловленность диагональных матриц Dp может быть довольно плохой. Для решения таких систем хорошо подходят прямые решатели разреженных систем. В качестве такого решателя был выбран PARDISO (PARallel Direct SOlver) из библиотеки Intel ® MKL [3]. Данный решатель является многопоточным, что позволяет реализовать гибридный подход и повысить производительность: на каждом используемом узле кластера можно решать задачи для одной подобласти, причем факторизацию СЛАУ, выполняющуюся при инициализации решателя, и решение серии систем с разными правыми частями можно выполнять в многопоточном режиме с использованием всех доступных процессорных ядер.
Для организации взаимодействия процессов на узлах кластера был использован интерфейс MPI, который широко используется при написании распределенных программ. При этом решатель был построен на основе master-slave архитектуры, с выделенным основным узлом (с рангом MPI, равным 0). На этом узле решателю
передается матрица, выполняется декомпозиция на подобласти, а также работает решатель GMRES. Остальные MPI процессы работают как ведомые, принимая от мастер-процесса матрицу системы на этапе инициализации и граничные переменные в процессе работы для решения задач в подобластях. Граничные для соседних подобластей переменные после решения отсылаются назад. Использование такого подхода позволяет на этапе итерационного решения ограничиться вызовами трех функций MPI на одну итерацию: MPI_scatter и MPI_gather для рассылки и сбора граничных переменных и MPI_Bcast для управления ведомыми процессами.
Представленный решатель был включен в состав библиотеки Krylov [5]. Интерфейс решателя, позволяющий вызывать его из процедур, написанных на различных языках программирования, представлен на рисунке.
void cschwarz_solver(MPI_Comm *pcomm, const int *n,
const double *a, const int *ia, const int *ja, double *u, const double *f, const double *eps, const int *ovl, const int *nmax, const int *nres, int *niter, double *time, int *err); void cschwarz_client(MPI_Comm *pcomm);
Интерфейс решателя Параметры решателя следующие:
pcomm — MPI коммуникатор, процессы которого участвуют в решении СЛАУ; процедура cschwarz_solver(...) должна вызываться в процессе с рангом 0, в остальных процессах необходимо вызывать cschwarz_client(...);
n — порядок СЛАУ;
a, ia, ja — матрица системы в формате CSR (Compressed Sparse Row format) [3];
u — вектор размера n, в который записывается полученное решение системы;
f — правая часть СЛАУ, должна иметь размер n;
eps — критерий останова итераций е (eps) для решателя GMRES, решение считается найденным, если
||g — [I — T]ub|| <e||g||;
ovl — параметр пересечения подобластей; после построения разбиения на подобласти без пересечений каждая подобласть расширяется на ovl фронтов влево и вправо, соответственно пересечение соседних подобластей состоит из 2 ■ ovl фронтов;
nmax — максимальное число итераций решателя GMRES;
nres — число итераций решателя GMRES до рестарта, позволяет ограничить длину последовательности qj базисных векторов подпространства Крылова; после выполнения каждых nres итераций, даже если не была достигнута необходимая
точность решения, текущее приближение пересчитывается по формулам
yj = eH-10 uj ^ u0 + [qo... qj-i]yj,
и новое значение ujj используется в качестве начального приближения для решателя;
niter — выходной параметр, количество совершенных решателем итераций; time — время, затраченное на решение системы;
err — код возврата решателя; ненулевое значение свидетельствует об ошибке.
3. Численные эксперименты
В последующих таблицах приняты следующие обозначения: N — порядок системы, Nnz — количество ненулевых элементов в исходной матрице, Nj — размерность пространства следов, tfac — время (в секундах) факторизации матриц Dp, ttot
— общее время решения и n — количество итераций. Эксперименты проводились на вычислительных узлах HP BL2х220с G6 кластера НКС-30Т [8], каждый из которых оборудован двумя четырехъядерными процессорами Intel Xeon Е5540 и 16-ю гигабайтами памяти. В целях увеличения объема памяти, доступной PARDISO для хранения факторов матриц Dp, запускалось по одному MPI-процессу на узел, а для ускорения факторизации были задействованы все имеющиеся на узлах ядра. Параметр пересечения подобластей при решении всех задач был установлен таким образом, чтобы подобласти пересекались максимум по 8 фронтам, за исключением второй модельной задачи, где пересечение составляет 4 фронта. Параметр е в критерии останова итераций был равен 10-7. Рестарты итерационного процесса не производились, то есть значение параметра nres было выбрано достаточно большим.
3.1. Модельная задача для диффузионно-конвективного уравнения
В качестве первой модельной задачи была выбрана задача Дирихле для линейного диффузионно-конвективного уравнения
Us + gf + §Zf + pfu + «Й + rf? = f (x.y.z), (x.y.z) e a
u|r = g(x,y,z),
в единичном кубе П = [0; 1]3 на равномерной кубической сетке. Дискретизация задачи проводилась с помощью семиточечной схемы экспоненциального типа [7]. Функции f и g выбирались соответственно точному решению u(x,y,z) = x2 + y2 + z2, а в качестве начального приближения была взята функция u0(x, y, z) = 0. В таблицах 1-3 приведены результаты решения данной модельной задачи при различных значениях p, q и r.
Максимальная погрешность решения первой модельной задачи в равномерной норме составила 7, 7 ■ 10-7.
Таблица 1
Решение первой модельной задачи при р = д = г = 0
Сетка N Количество узлов
2 4 8 16
6 со 262144 1810432 N п ^!ае 6104 12 3,31 4,65 16196 16 0,94 1,74 34688 21 0,46 1.18 71786 32 0,30 1,26
1283 2097152 14581760 Nb п tfac 24536 17 219 245 65508 23 33,3 44,5 139183 29 8,36 15,6 285485 40 2,56 9,25
Таблица 2
Решение первой модельной задачи при р = д = г = 16
Сетка N Количество узлов
2 4 8 16
6 со 262144 1810432 N п tfac 6104 10 3,29 4,44 16196 11 0,95 1,67 34688 14 0,50 1,11 71786 24 0,30 1,08
1283 2097152 14581760 N п tfac 24536 14 233 255 65508 16 34,7 43,9 139183 20 8,73 14,8 285485 29 2,60 8,23
Таблица 3
Решение первой модельной задачи при р = д = г = -16
Сетка N Количество узлов
2 4 8 16
6 со 262144 1810432 N п tfac 6104 10 3,33 4,48 16196 11 1,00 1,97 34688 15 0,47 1,14 71786 25 0,44 1,33
1283 2097152 14581760 N п tfac 24536 14 226 248 65508 17 33,9 43,5 139183 21 8,67 15,2 285485 30 3,01 8,77
3.2. Модельная задача с волноводом
Вторым примером является решение трехмерной краевой задачи для уравнения Гельмгольца
V х | — V х E j — E = 0, fcg = wJ£0^0,
\Vr )
дискретизация которой выполнялась с помощью конечно-элементных аппроксимаций с базисными функциями Неделека второго порядка [9, 10]. Расчетной областью является волновод с линейными размерами a = 72, b = 34, c = 200 мм. Параметры
задачи: = 1, £r = 1 — 0, li, частота w = 6п • 109 Гц. Граничные условия: на грани
z = 200 мм задавалась касательная компонента поля
Е0 х П = ey sin(nx/a) х n,
а на остальных гранях — условие металлической стенки (E0 = 0). Аналитическое решение в данном случае имеет вид
^ / nx \ sin yz
E = ey sin — -------------.
V a / sin yc
В табл. 4 приводятся результаты решения указанной задачи при различных размерностях сетки.
Таблица 4
Решение второй модельной задачи
Сетка N Nnz Количество узлов
2 4 8 16
8 x 4 x 20 21664 423392 Nb n tfac ttot 1540 9 0,19 0,51 4678 13 0,17 0,47 10708 21 0,15 0,67 23052 258 0,36 6,19
15 x 7 x 40 149874 3117779 Nb n tfac ttot 5210 13 2,56 5,88 16646 18 1,26 4,00 37670 25 0,63 3,12 78446 38 0,40 3,15
29 x 14 x 80 1196026 25795767 Nb n tfac ttot 20562 18 108 173 67084 26 39,9 85,6 151032 34 14,9 43,3 318108 49 5,29 32,94
Наличие точного решения позволяет проверить сходимость итерационного процесса. Данные тесты продемонстрировали ошибку порядка 0(Л,2), что согласуется с теорией.
3.3. Примеры решения практических задач
В заключение данного пункта мы приведем результаты решения трех СЛАУ (обозначаемых ниже в табл. 5 как Г1, Г2, Г3), возникающих из актуальных практических краевых задач, описываемых трехмерными уравнениями гидродинамики, которые аппроксимируются с помощью неявных конечно-объемных схем на сложных неструктурированных сетках.
Таблица 5
Решение практических задач гидродинамики
Задача N Количество узлов
2 4 8 16
Г1 4875172 26455278 N п ^/ас 180487 40 101 166 439384 58 64 117 895041 85 49,8 119 1831917 162 19,7 201
Г2 1726272 11984539 N п І/ас 6912 52 54,4 79,0 18533 109 26,0 116 42476 207 9,55 83,0 90000 334 3,26 64,0
Г3 475000 13672497 N п /ас Ііоі 2500 24 6,12 10,1 7500 25 2,62 5,45 17500 36 1,08 3,31 37110 53 0,53 2,88
Следует отметить, что приведенные результаты параллельных вычислений дают значительный выигрыш в сравнении с однопроцессорными расчетами, в которых некоторые из рассмотренного типа СЛАУ с помощью «стандартных» предобуслов-ленных крыловских алгоритмов решаются часами.
Заключение
Рассмотренные результаты позволяют сделать следующие выводы.
• С увеличением количества подобластей число итераций увеличивается, но (за одним исключением) не более чем линейно. В этом плане улучшение сходимости можно ожидать за счет включения других типов граничных условий на смежных границах, вместо используемого условия Дирихле.
• Применение прямого решателя РАИДКО для подобластей является достаточно экономичным. В частности, с увеличением числа МР1-процессов относительные временные расходы на факторизацию вспомогательных подсистем значительно уменьшаются.
• Размерности СЛАУ в подпространствах следов (N5) значительно меньше исходных порядков N. Это позволяет экономично реализовать с помощью GMR.ES итерации по подобластям.
• Описанная программная реализация позволяет проводить оптимизацию алгоритмов за счет выбора различных счетных параметров (величина перехлестов, варьирования числа ядер на подобласти, количества MPI-процессов на вычислительный узел и т.д.).
• При увеличении числа подобластей и используемых процессоров до нескольких тысяч целесообразен, по-видимому, переход от используемой 1-D декомпозиции к 2-D или даже 3-D декомпозиции.
Работа выполнена при поддержке грантов РФФИ № 11-01-00205 и Президиума РАН № 2.5.
Литература
1. Karypis, G. A Fast and Highly Quality Multilevel Scheme for Partitioning Irregular Graphs / G. Karypis, V. Kumar // SIAM Journal on Scientific Computing. - 1999. -Vol. 20, № 1. - P. 359-392.
2. Saad, Y. Iterative Methods for Sparse Linear Systems, Second Edition. / Y. Saad. -SIAM, 2003.
3. Intel (R) Math Kernel Library from Intel: сайт. URL: http://software.intel.com/ en-us/articles/intel-mkl/ (дата обращения: 03.11.2012).
4. Eigen: сайт. URL: http://eigen.tuxfamily.org/index.php?title=Main_Page
(дата обращения: 03.11.2012).
5. Ильин, В.П. Krylov: библиотека алгоритмов и программ для решения СЛАУ /
B.П. Ильин, Д.С. Бутюгин, Е.А. Ицкович и др. // Современные проблемы математического моделирования. Математическое моделирование, численные методы и комплексы программ. Сборник трудов Всероссийских научных молодежных школ. - Ростов-на-Дону: Изд-во Южного федерального университета, 2009. -
C. 110-128.
6. Кормен, Т. Алгоритмы: построение и анализ / Т. Кормен, Ч. Лейзерсон, Р. Ри-вест. - М: МЦНМО: БИНОМ. Лаборатория знаний, 2004.
7. Ильин, В.П. Методы и технологии конечных элементов. / В.П. Ильин. - Новосибирск: Изд-во ИВМиМГ СО РАН, 2007.
8. Кластер НКС-30Т: сайт. URL: http://www2.sscc.ru/HKC-30T/HKC-30T.htm (дата обращения: 03.11.2012).
9. Monk, P. Finite Element Methods for Maxwell’s Equations. / P. Monk. - Oxford University Press, 2003.
10. Ingelstrom, P. A new set of H(curl)-conforming hierarchical basis functions for tetrahedral meshes / P. Ingelstrom // IEEE Transactions on Microwave Theory and Techniques. - 2006. - Vol. 54, № 1. - P. 160-114.
Дмитрий Сергеевич Бутюгин, младший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН, аспирант, Новосибирский Государственный Университет, dm.butyugin@gmail.com.
Валерий Павлович Ильин, доктор физико-математических наук, профессор, главный научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН, ilin@sscc.ru.
Данил Валерьевич Перевозкин, младший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН, dperevozkin@mail.ru.
PARALLEL METHODS FOR SLAE SOLUTION ON THE SYSTEMS WITH DISTRIBUTED MEMORY IN KRYLOV LIBRARY
D.S. Butyugin, Institute of Computational Mathematics and Mathematical Geophysics SB RAS (Novosibirsk, Russian Federation),
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 presents an approach to creation of black-box parallel iterative solver, which is used in Krylov library for solving systems of linear algebraic equations (SLAEs) with large sparse matrices in CSR format arising from discretization of multidimensional boundary value problems. A variant of one-dimensional algebraic decomposition method is offered. The algorithm is based on breadth-first search of SLAE’s adjacency graph that allows to reduce the matrix to block-tridiagonal form. The algebraic solver is based on additive Schwarz method which naturally suits distributed memory computer systems. The generalized minimal residual method is used to solve the SLAEs arising from relations on subdomains’ boundaries. Auxiliary subdomain systems are solved with Intel MKL’s multithreaded direct solver PARDISO. Implemented algorithms were tested on the numerical solution of the series of computational mathematics problems, such as problems of hydrodynamics, diffusion-convection equations, problems of electromagnetism and others. Adduced numerical experiments results show the effectiveness of the presented algorithms for multiprocessor computational systems with distributed memory..
Keywords: iterative algorithms, domain decomposition methods, parallel computing, algebraic systems, sparse matrices, numerical experiments, additive Schwarz method.
References
1. Karypis G., Kumar V. A Fast and Highly Quality Multilevel Scheme for Partitioning Irregular Graphs. SIAM Journal on Scientific Computing. 1999. Vol. 20, № 1. P. 359392.
2. Saad Y. Iterative Methods for Sparse Linear Systems, Second Edition. Society for Industrial and Applied Mathematics, 2003.
3. Intel (R) Math Kernel Library from Intel.
URL: http://software.intel.com/en-us/articles/intel-mkl/
4. Eigen. URL: http://eigen.tuxfamily.org/index.php?title=Main_Page
5. Il’in V.P., Butyugin D.S., Itskovich E.A et al. Krylov: biblioteka algoritmov
i programm dlya resheniya SLAU [Krylov: a Library of Algorithms and
Programs for Solving SLAE]. Trudy Vserossiyskikh nauchnykh molodezhnykh shkol ”Sovremennye problemy matematicheskogo modelirovaniya. Matematicheskoye modelirovaniye, chislennye metody i kompleksy programm” [Proceedings of All-Russian Youth Scientific Schools ”Modern Problems of Mathematical Modeling. Mathematical Modeling, Numerical Methods and Program Complexes]. Rostov-on-Don, SFedU Publ., 2009, P. 110-128.
6. Cormen T., Leiserson C., Rivest R. Introduction to Algorithms. MIT Press, 2001.
7. Il’in V.P. Metody i tekhnologii konechnykh elementov [Methods and Technologies of Finite Elements]. Novosibirsk, ICM&MG SBRAS Publ., 2007.
8. NKS-30T Cluster URL: http://www2.sscc.ru/HKC-30T/HKC-30T.htm
9. Monk P. Finite Element Methods for Maxwell’s Equations. Oxford University Press, 2003.
10. Ingelstrom P. A New Set of H(curl)-conforming Hierarchical Basis Functions for Tetrahedral Meshes IEEE Transactions on Microwave Theory and Techniques. 2006. Vol. 54, № 1. P. 160-114.
Поступила в редакцию 5 ноября 2012 г.
