CC BY
60
УДК 519.612
ФУНКЦИОНАЛЬНОСТЬ И ТЕХНОЛОГИИ АЛГЕБРАИЧЕСКИХ РЕШАТЕЛЕЙ В БИБЛИОТЕКЕ KRYLOV1
Д.С. Бутюгин, Я.Л. Гурьева, В.П. Ильин, Д.В. Перевозкин, А.В. Петухов, И.Н. Скопин
Описываются функциональные возможности и особенности программной реализации библиотеки параллельных алгоритмов Krylov, ориентированной на решение больших систем линейных алгебраических уравнений с разреженными симметричными и несимметричными матрицами (положительно определенными и знаконеопределенными), получаемых при сеточных аппроксимациях многомерных краевых задач для систем дифференциальных уравнений на неструктурированных сетках. Библиотека включает двухуровневые итерационные методы в подпространствах Крылова, предобуславливание которых осуществляется на основе сбалансированной декомпозиции расчетной области с различными размерами пересечений подобластей и краевых условий сопряжения на смежных границах. Программные реализации выполнены на типовых сжатых разреженных форматах матричных данных. Приводятся результаты численных экспериментов с демонстрацией эффективности распараллеливания для характерных плохо обусловленных задач. Ключевые слова: предобусловленные итерационные алгоритмы, подпространства Крылова, методы декомпозиции областей, разреженные алгебраические системы, численные эксперименты.
Введение
Настоящая работа содержит описание функционального наполнения и технологических подходов библиотеки параллельных итерационных алгоритмов Krylov (см. предварительные публикации в [1, 2]), ориентированной на решение больших систем линейных алгебраических уравнений (СЛАУ) с разреженными матрицами, возникающими при конечнообъемных или конечно-элементных (МКО или МКЭ [3, 4]) аппроксимациях многомерных краевых задач для систем дифференциальных уравнений на неструктурированных сетках, на многопроцессорных вычислительных системах (МВС с количеством ядер в десятки и сотни тысяч). В данном случае имеется в виду отсутствие программных ограничений на проведение таких масштабных вычислительных экспериментов, а также достаточно высокая эффективность применяемых алгоритмов.
Основой применяемых в библиотеке Krylov вычислительных подходов являются двухуровневые итерационные процессы в подпространствах Крылова, предобуславливаемые с помощью аддитивного метода Шварца и декомпозиции расчетной области с пересечениями подобластей и разными типами краевых условий на смежных внутренних границах, см. [4-7] и цитируемые там работы.
Внешний итерационный процесс осуществляется распределенным по вычислительным узлам образом методом FGMRES [7] с динамическими (в общем случае) предобуславливате-лями, которые включают решение вспомогательных подсистем в расширенных подобластях с помощью или прямого решателя PARDISO из библиотеки Intel MKL [8], или авторскими
1 Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2013»
версиями алгоритмов BiCGStab и GMRES [7, 9], предобусловленных с помощью покомпонентной или эффективной «мелкоблочной» модификациями Айзенштата неполной факторизации. В последнем случае имеются в виду СЛАУ, полученные, например, при сеточных аппроксимациях систем дифференциальных уравнений, когда одному узлу сетки соответствует несколько неизвестных функций (в задачах тепломассопереноса это могут быть три компоненты вектора скорости, давление и температура). Отличительной чертой соответствующих матриц является то, что они представимы в блочном виде с диагональными матрицами малого порядка, который не зависит от числа узлов сетки и, следовательно, от порядка СЛАУ.
Коммуникации между процессорами на внешних итерациях выполняются средствами системы MPI, а «внутренние» вычисления при решении СЛАУ в подобластях реализуются параллельными потоками над общей памятью с помощью OpenMP. Библиотека Krylov включает также итерационные решатели из написанной на языке CUDA библиотеки Cusp [10], что позволяет решать на гетерогенном кластере внутренние подсистемы на
Особенностью программных реализаций является то, что решаемые СЛАУ представляются в стандартном сжатом строчном формате CSR [8], который является практически безальтернативным способом построения универсальных алгебраических решателей, но представляет значительные трудности для эффективного распараллеливания алгоритмов. В частности, возникает нетривиальная задача формирования CSR-форматов для вспомогательных расширенных СЛАУ в подобластях.
Библиотека Krylov является функциональным аналогом таких пакетов распределенных итерационных решателей, как например PETSc [11], Hypre [12], HIPS и pARMS [13] и др., однако в Krylov реализован также ряд оригинальных методов, обсуждаемых в работе.
В разделе 1 мы описываем функциональное наполнение библиотеки Krylov, в разделе 2 излагаются особенности программных реализаций алгебраических решателей, последний же раздел посвящен численным экспериментам для характерных практических задач. В заключении подводятся итоги работы и делаются выводы о применимости изложенных методов.
1. Функциональное наполнение библиотеки Krylov
Рассматривается решение вещественной СЛАУ
матрица которой, возможно, предварительно разбита на блочные строки А = {Ар є ^мр,и, р = 1,..., р}, N + ... + N = Ж, с приблизительно равным числом строк N ^ 1 в каждой. Соответствующим образом также разбиваются векторы искомого решения и пра-
GPU.
Au = f, A = {aij} є Rn,n,
(І)
вой части и = {ир}, / = {/р}, так что система уравнений (1) может быть записана в форме совокупности Р подсистем меньших порядков:
р
(2)
q = 1
q=p
где Ар,д є 'Я,Мр’Мі — прямоугольные матричные блоки, получаемые при разбиении каждой блочной строки (и всей матрицы А) на блочные столбцы.
Будем считать, что СЛАУ (1) представляет собой систему сеточных уравнений, так
что каждая компонента векторов u, f соответствует узлу сетки, общее число которых в
P
расчетной сеточной области Qh = У Qp равно N. При этом блочное разбиение матриц и
р=1
векторов соответствует разбиению (декомпозиции) Qh на P сеточных непересекающихся подобластей Qp, в каждой из которых находится Np узлов.
Описанная декомпозиция Qh не использует узлов-разделителей, т.е. условные границы подобластей не проходят через узлы сетки. Формальности ради можно считать, что внешность расчетной области есть неограниченная подобласть Qo с числом узлов No = 0. Сеточное уравнение для i-го узла сетки может быть записано в виде
где через Шг обозначается сеточный шаблон г-го узла, т.е. совокупность номеров всех узлов, участвующих в г-м уравнении.
Предполагается, что сеточные узлы и соответствующие переменные пронумерованы следующим образом: сначала идут подряд все узлы 1-й подобласти Оі (неважно в каком внутреннем порядке), затем — узлы второй подобласти и т.д. Отметим, что взаимнооднозначного соответствия алгебраической и геометрической (сеточной) интерпретации структуры СЛАУ может и не существовать. Типичный пример — одному узлу сетки может соответствовать т > 1 переменных в алгебраической системе. Если для каждого узла такие «кратные» переменные нумеруются подряд, то структура СЛАУ приобретает «мелкоблочный» (т ^ N) вид (в (3) иг и / означают не числа, а подвекторы порядка т соответственно, аг^ Є ҐКт,т) и на таких случаях мы остановимся в последующем особо.
Пусть матрица (1) задана в сжатом разреженном СБИ-формате [8] с указанием числа строк Жр в каждой из Р подсистем (2), в соответствии с разбиением матрицы А на блочные строки Ар. Естественно, предполагается, что строки исходной матрицы пронумерованы
р-1 р
подряд от 1 до N, так, что номера строк каждого блока Ар меняются от 1 + ^ N до ^ N.
На основе блочного представления СЛАУ (2) стандартным образом строится аддитивный метод Шварца, на алгебраическом языке представляющий блочный итерационный метод Якоби. Однако известно, что скорость сходимости итераций возрастает, если декомпозицию расчетной области сделать с пересечением областей.
В силу этого мы используем определение расширенной сеточной подобласти Ор Э Пр, имеющей пересечения с соседними подобластями, величину которых мы будем определять в терминах количества околограничных сеточных слоев, или фронтов. Обозначим через Г0 £ Ор множество внутренних околограничных узлов из Ор, т.е. таких узлов Рг £ Ор, у которых хотя бы один из соседей не лежит в &р(Р] £ &р, ] £ Шг, ] = г)). Обозначим далее через Гр множество узлов, соседних с узлами из Гр, но не принадлежащих Ор, через Гр — множество узлов, соседних с узлами из Гр, но не принадлежащих объединению Гр и Ор, через Гр — множество узлов, соседних с Гр, но не принадлежащих Гр и Гр и Ор, и т.д. Соответственно эти множества назовем первым внешним слоем (фронтом) узлов, вторым слоем, третьим и т.д. Получаемое объединение узлов
(З)
(4)
будем называть расширенной p-й сеточной подобластью, а целую величину А определяем как величину расширения, или пересечения (в терминах количества сеточных слоев). Случай А = 0 фактически означает декомпозицию области Qh на подобласти без пересечений (Q0 = Qp).
На формальном алгебраическом языке каждой расширенной подобласти можно сопоставить подсистему уравнений
P
(Ap,p + eDp)Up = fp 'У ' Ap,qUq + ODpUp = f p, (5)
q = i q=p
где Up, fp, rp £ RNp, а Dp — диагональная матрица, определяемая соотношением
P
Dpe = X! Ap,q e, e = (1,..., 1)T eRiVp. (6)
q = 1 q=p
Здесь в £ [0,1] — некоторый итерационный параметр, который при в = 0 соответствует итерационному краевому условию Дирихле на смежных границах подобластей, при в = 1 — условию Неймана, а при 0 < в < 1 — условию 3-го рода, или Робена.
Переходя теперь к полному вектору и и вводя матрицы
Bp = Ap,p + вDp, (7)
из (5) получаем итерационный аддитивный процесс Шварца в виде
ип = un-1 + B-1(f - Aun-1) = un-1 + B-1rn, (8)
где B — предобуславливающая матрица, определяемая следующим образом:
P
B-1 = B-S = Е Bp-1, Bp-1 = WpTBp-1Wp. (9)
p=i
Здесь Wp : Rn ^ RNp есть матрица сужения полного вектора в подвектор из Qp, а Wp^ — транспонированная матрица продолжения (расширения вектора из Qp в Q).
Выполнение каждой итерации в (8) требует параллельного решения внутренних СЛАУ в подобластях Qp, что может производиться или прямым методом, или итерационным алгоритмом крыловского типа. В последнем случае фактически реализуются переменные (динамические) предобуславливатели Bn, при этом реально вместо (8) используется универсальный «гибкий» метод обобщенных минимальных невязок FGMRES [7], оптимизирующий невязку в подпространствах Крылова.
2. Особенности параллельной реализации алгоритмов
Описанные выше параллельные алгоритмы реализованы в форме библиотеки Krylov, которая ориентирована на решение сверхбольших СЛАУ на системах с распределенной памятью, включающих в себя большое количество вычислительных узлов. В связи с этим реализация итерационных решателей СЛАУ и предобуславливателей имеет определенные особенности. В частности, требуется такая организация и структура методов, которая хорошо отображается на архитектуры имеющихся вычислительных систем.
2.1. Общая организация двухуровневых итераций
Решатели в библиотеке Krylov построены на основе широко распространенного стандарта интерфейса обмена данными MPI (Message Passing Interface). Базой для итерационного решения систем уравнений в библиотеке Krylov является метод FGMRES, распределенный по различным MPI-процессам. Данный метод совместно с различными предобуславливате-лями типа Шварца используется для решения заданной распределенной СЛАУ.
Псевдокод алгоритма представлен на рисунке.
г0 4-f-Auo, Р = ||г0||, qo «-r0/||r0||, £ «- (1,0,... ,0)т For n 6 {0,1,... } while ||rn|| > e/3
%Т1 ^ Bn (]n
Qn+l ^ Azn
For fee [0,...,n]
Hk,n 4 {.Qn+hQk) Qn+l 4 Qn+l Hfc^nQk
EndFor
Hn+l,n 4 | |Qn+l | 15 Qn+l 4 9n+l/-^n+l,n
For к E {0,..., n — 1}
1 1 2 Ск Sk 1 £ 3 1
-^fc+l,n л ~8к Ск Hk+\,n
EndFor
Cn ^ \Hn,n\/-y/|Hn,n^ “Ь |-Ип+1,пР
&n ^ СпНп+1,п/Hn,n
1 *5 | 2. 1 а? dr I 1 i
1 гН + *5 | л 1 dr i«o 1 1 £n+l
Нп,п ^ Сп-Нг»,П "Ь &пНп+ 1,71) Нп+1,п ^ 0
I |гтг+111 <- /З|^гг-Ы |
EndFor Уп <- /ЗН-^
Хп <- ж0 + [г0 ... гп-1 \уп.
Метод FGMRES
Анализ вычислительной схемы алгоритма показывает, что основными в методе FGMRES являются следующие операции:
• векторно-векторные операции (вычисление скалярных произведений, норм, и т.п.);
• матрично-векторные операции (умножение матрицы на вектор); в первых двух пунктах операции выполняются в полном N-мерном пространстве, где N — размерность решаемой СЛАУ;
• применение предобуславливателя (метод Шварца);
• ^Я-разложение матрицы Хессенберга, формируемой в процессе построения базисных векторов.
Все остальные операции требуют 0(1) времени и не вносят существенного вклада в производительность решателя. Более того, вклад ^Я-разложения матрицы Хессенберга во время работы решателей также можно игнорировать в предположении, что порядок исходной системы много больше количества итераций. Такое разложение также не требует
большого количества памяти, поэтому с целью уменьшения коммуникационных затрат оно выполняется в каждом МР1-процессе.
Векторно-векторные операции в итерационных решателях распараллеливаются естественным образом за счет разбиения векторов, порождаемого декомпозицией матрицы на расширенные блочные строки. Единственная особенность реализации этих операций заключается в том, что для вычисления скалярного произведения и нормы необходима будет точка синхронизации и обмен данными, однако, к примеру, в интерфейсе МР1 имеется требуемая функция МР1 АПге^се, реализация которой оптимизируется поставщиком библиотеки МР1.
Для проведения матрично-векторных операций решатель строит для каждой подобласти списки узлов, являющихся граничными для соседних подобластей. В дальнейшем, при умножении матрицы на вектор или при применении итерации Шварца к вектору неизвестных, решатель использует полученную информацию для пересылки частей вектора между процессами. Ключевым моментом в таком подходе является уменьшение объема пересылаемых данных — необходимо пересылать в соседние подобласти только граничные элементы вектора вместо того, чтобы пересылать вектор целиком. В случае, если разбиение на подобласти было построено подходящим способом, количество граничных вершин Nг будет существенно меньше размера подобласти (например, для двумерных задач N будет пропорционально квадратному корню из общего числа вершин в подобласти, а для трехмерных N ~ N2/3).
Применение предобуславливателей, основанных на методе Шварца, требует на каждой итерации внешнего алгоритма решения подзадач в подобластях. Достоинством метода аддитивного Шварца является то, что решение в подобластях может проводиться независимо друг от друга и не требует коммуникаций. Действительно, как уже отмечалось выше, итерацию Шварца можно представить в виде (8), где предобуславливатель В составлен из блоков (7). Отсюда видно, что применение предобуславливателя В-1 в подобласти с номером р не требует знания переменных, не принадлежащих ей.
2.2. Реализация внутренних итераций в подобластях
Для решения задач Ар,рур = Ьр в подобластях можно использовать какой-нибудь прямой решатель для разреженных систем, например решатель PARDISO из библиотеки 1гЛе1® МКЬ. Однако время, требуемое PARDISO для разложения матриц, в общем случае растет практически как О^3/Р3), поэтому такой метод подходит только для решения не очень больших систем на большом числе процессоров. Однако действие обратных несимметричных матриц А—р можно вычислять приближенно при помощи предобусловленных методов в подпространствах Крылова, таких как GMRES и BiCGStab. В качестве предобуславливателя предлагается метод USOR в модификации Айзенштата, в том числе его мелкоблочный вариант. При этом пользователю предоставляется возможность выбора решателя для подобластей — либо PARDISO из 1гЛе1® МКЬ, либо одного из перечисленных выше итерационных решателей.
Рассмотрим более подробно предобуславливатель USOR в модификации Айзенштата. Пусть А = О — Ь — и, а В — предобуславливающая матрица, записанная в форме
В = (С — Ь)С-1(С — и ) = С — Ь — и + ЬС-1 и, (10)
где G = block-diag{ Gk}, Gk G 'Rm,m, k = = N/m есть блочно-диагональная
невырожденная матрица с диагональными блоками одинакового порядка m, причем порядок СЛАУ кратен m, а M — это блочный порядок матрицы A. Матрица G — некоторая пока не конкретизируемая невырожденная матрица, для которой вычислимо разложение G = LqUg с треугольными невырожденными множителями Lg и Ug.
С помощью обозначений
Важно отметить, что умножение А на некоторый вектор у можно представить в удобной для вычисления форме
где ш — релаксационный числовой параметр. Очевидно, что такая матрица G сохраняет блочно-диагональную структуру В. В этом случае треугольное разложение G сводится к соответствующему разложению ее диагональных блоков:
а реализация умножения вектора на матрицу А упрощается за счет блочного представления матриц
Точнее, при использовании блочного представления векторов у = {ук}, и = {и к}, к = 1, ...,М, решения участвующих в (14) треугольных систем
D = L-1DU-1, L = L-1LU-1, U = LG1UU-1
(ll)
предобусловленная матрица системы записывается в форме
U = (I - L)-1L^1AU-1(I - U)-1 =
= (I - L)-1 + (I - U)-1 + (I - L)-1(D - 2I)(I - U)-1,
(12)
где I есть единичная матрица.
Отсюда предобусловленную СЛАУ можно представить как
Au = f = Lg1(I - L) 1f, и = (I - U)UGu.
(13)
Av = (I - L) l[v + (D - 2I)w] + w, w = (I - U) 1v,
(14)
составляющей суть модификации Айзенштата.
Конкретизируем теперь выбор матрицы С. Мы используем
(15)
G = block-diag{ Gk = LG,k UG,k },
(16)
L = {Lki = LG1kLk,iU-1}, U = {Uk,i = l~g\UkiU-1}, k,l = 1,...,M.
(17)
(I - U)w = v, (I - L)y = z = v + (D - 2I)w
(1В)
осуществляются фактически по следующим рекуррентным формулам:
м
wm = vm, wk = vk + ^ Uk,lvi, k = M - 1,..., 1, l=k+1
k-1
(19)
Отметим, что используемые модификации неполной факторизации являются экономичными, но плохо распараллеливаемыми, вследствие необходимости решения вспомогательных СЛАУ с треугольными матрицами. Однако имеются различные подходы к распараллеливанию треугольных решателей, например метод планирования вычислений по уровням, см. [7]. Кроме того, умножение на «мелкоблочные» матрицы Uk,i и Lk,i, если они вычислены заранее до итераций, для каждого фиксированного k позволяют обеспечить существенное ускорение средствами OpenMP.
3. Результаты численных экспериментов
Мы приведем результаты некоторых численных экспериментов по решению ряда практических задач с помощью библиотеки Krylov, выполненных на кластере ИВМиМГ СО РАН [14]. В данных расчетах внешние итерации проводились с помощью метода FGMRES с критерием окончания итераций по условию на евклидовую норму невязки
1И|2 < е||/||2, е = 10-7. (20)
Начальные приближения для искомых векторов всегда выбирались и0 = 0.
Вспомогательные СЛАУ в подобластях решались либо с помощью прямого решателя PARDISO, причем наиболее трудоемкий этап — LU-разложение матрицы — выполнялся один раз до итераций, либо с помощью итерационного метода BiCGStab с предобуславли-вателем Айзенштата. Расчеты проводились на различном количестве P вычислительных узлов, с формированием такого же числа подобластей и соответствующих им MPI-процессов. В нижеследующих таблицах N и NZ обозначают порядок решаемых СЛАУ, а также количество ненулевых элементов матрицы, которые характеризуют трудоемкость задачи.
Рассмотрим результаты экспериментов по решению трех несимметричных СЛАУ, возникающих из сеточных аппроксимаций практических задач гидро-газодинамики. Характеристики соответствующих плохо обусловленных матриц даны в табл. 1.
Таблица І
Характеристики «гидро-газодинамических» матриц
Наименование 6 - 0 1 NZ • 10-6
Г2 1,73 12,0
Г3 0,4В 13,7
Г4 9,3В 50,4
В табл. 2 — 4 приведены результаты расчетов для данных СЛАУ с использованием PARDISO в качестве решателя в подобластях. В каждой из клеток таблиц сверху вниз представлено общее время счета t в секундах и количество внешних итераций ne.
Здесь и далее используются следующие обозначения: Nnod — число используемых вычислительных узлов; Nmpi — количество MPI-процессов на одном узле; P = Nnod • Nmpi — общее число MPI-процессов, равное количеству подобластей; А — параметр пересечения в декомпозиции областей (количество расширений сеточной подобласти); Nth — количество формируемых вычислительных потоков в одном MPI-процессе; Ti — время (в секундах) решения СЛАУ без распараллеливания, т.е. при Nmpi = Nth = 1.
Таблица 2
Решение СЛАУ Г2 с помощью РАИБШО (Т\ = 87,1)
P \ А 0 1 2 3 4 5
5 60,6 235 39,9 116 33,9 78 30,9 59 28,5 46 27,4 37
10 58,4 486 32,9 244 24,6 167 21,0 124 19,1 99 17,4 79
20 118,5 934 47 479 28,6 330 22,6 246 18,3 194 16,0 157
30 230,5 1352 74,5 690 53,9 473 33,4 355 27,0 280 20,8 228
Таблица 3
Решение СЛАУ Г3 с помощью РАИБШО, Щь, = 4 (Т\ = 10,1)
P \ А 0 1 2 3 4 5
5 9,45 82 7,46 41 7,23 29 6,80 22 6,75 18 6,52 15
10 5,64 114 4,69 59 4,32 41 4,21 32 4,27 27 4,31 23
20 6,03 164 4,09 84 3,72 58 3,76 44 3,77 37 3,59 32
30 6,83 183 4,29 94 3,86 67 3,63 52 3,78 43 3,68 38
Таблица 4
Решение СЛАУ Г4 с помощью РАИБШО (Т\ = 399, 4)
P \ А 0 1 2 3 4 5
5 280,4 235 199,4 119 153,9 85 146,5 67 140,6 55 143,5 47
10 127,9 311 73,6 161 61,3 116 56,3 92 51,9 76 49,6 66
20 147,4 331 71,2 175 51,9 127 47,0 101 41,3 84 40,6 72
30 158,3 355 76,9 186 59,2 135 52,7 109 49,5 95 48,4 80
Как видно из этих результатов, использование пересечений с ростом А дает стабильное уменьшение числа внешних итераций, а время счета уменьшается в два и в большее число раз (максимальный коэффициент ускорения — свыше 10). С ростом Р время вычислений сначала уменьшается, но после Р > 20 начинает увеличиваться, так как растет число внешних итераций. Для исправления ситуации, очевидно, надо использовать ускорение типа грубосеточной коррекции [15], которое в данных экспериментах не применялось.
Второй набор СЛАУ для численных экспериментов был представлен в блочно-строчном распределенном формате СБИ, а заданное количество блочных строк указано в табл. 5. Там же указано, что первые две СЛАУ данного набора имеют «мелкоблочную» структуру с размером блоков, равным т = 5. В таблице также приведено количество подобластей Р, на которые были разбиты соответствующие СЛАУ. Разбиение СЛАУ осуществлено без пересечений, так что, фактически, в данном наборе тестов А = 0.
Таблица б
Характеристики блочных распределенных СЛАУ
Наименование N•lO-6 NZ • lO-6 P m
Г5 9,5 330 48 5
Г6 3,7 126,8 36 5
Г7 6,5 44,4 20 1
Г8 0,35 1,74 20 1
В табл. 6 приведены результаты экспериментов для четырех СЛАУ из второго набора тестируемых задач. В клетках таблицы указано общее время счета Ь (слева), а также число внешних итераций пе.
Для внутренних СЛАУ при использовании итерационных решателей фактически использовались ограничения числа итераций пгтах. Для прямого решателя РАКБШО в скобках указано число запускаемых им потоков при факторизации и решении СЛАУ. В качестве предобуславливателя для внутренних итерационных решателей использовался «мелкоблочный» предобуславливатель ИБОИ в модификации Айзенштата. Для предобуславливателя ИБОИ представлены результаты без использования средств распараллеливания ОрепМР.
Таблица 6
Сравнение прямых и итерационных решателей для СЛАУ в подобластях
Метод Г5 Г6 Г7 Г8
PARDISO (Nth = S) 58,9 38 18,9 21 121,6 351 3,8 189
PARDISO (Nth = 2) 149,6 38 42,4 21 178,8 351 4,7 189
BBiCGStab (n%max = lO) 49,3 39 13,6 21 498,5 368 14,6 252
BBiCGStab (nlmax = 2O) 70,5 38 19,7 21 843,9 360 24,9 206
Полученные результаты позволяют сделать следующие выводы:
• на рассматриваемых СЛАУ с m = l прямой решатель PARDISO имеет существенное преимущество перед итерационными, причем увеличение в нем количества используемых потоков Nth от 2 до 8 дает коэффициент ускорения от 1,5 до 3;
• прямой решатель PARDISO проигрывает «мелкоблочному» итерационному алгоритму для СЛАУ с m = 5, хотя текущая реализация даже не использует средства OpenMP для распараллеливания на общей памяти;
• расчеты подтверждают ожидаемый факт, что при использовании итерационных внутренних решателей для СЛАУ в подобластях нет смысла добиваться высокой точности на различных внешних итерациях; более конкретно, при увеличении количества внут-
ренних итераций n%max с 10 до 20 число внешних итераций несколько падает, но каждая из них делается «дороже» и общее время решения увеличивается.
Заключение
В статье рассмотрен параллельный двухуровневый итерационный решатель, входящий в состав библиотеки Krylov, и приведены ключевые моменты его реализации: организация двухуровневых итераций с использованием методов Шварца и FGMRES, решение СЛАУ в подобластях итерационными и прямыми методами. Произведены эксперименты, показывающие применимость изложенных подходов к решению характерных задач. В частности, продемонстрирована возможность использования итерационных методов с невысокой точностью для решения СЛАУ в подобластях, а также показано преимущество таких методов в случае явного учета блочной структуры матриц.
Работа поддержана грантом РФФИ N 11-01-00205, а также грантами Президиума РАН N 15.9-4 и ОМН РАН N 1.3.3-4.
Литература
1. Krylov: библиотека алгоритмов и программ для решения СЛАУ / Д.С. Бутюгин,
В.П. Ильин, Е.А Ицкович и др. // Современные проблемы математического моделирования. Математическое моделирование, численные методы и комплексы программ. Сборник трудов Всероссийских научных молодежных школ. Ростов-на-Дону: Изд-во Южного федерального университета, 2009. — С. 110-128.
2. Бутюгин, Д.С. Методы параллельного решения СЛАУ на системах с распределенной памятью в библиотеке Krylov / Д.С. Бутюгин, В.П. Ильин, Д.В. Перевозкин // Вестник ЮУрГУ. Серия «Вычислительная математика и информатика». — 2012. — №. 47(306).
— С. 5-19.
3. Ильин, В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений / В.П. Ильин — Новосибирск: Изд-во ИВМиМГ СО РАН, 2001. — 345 с.
4. Ильин, В.П. Методы и технологии конечных элементов / В.П. Ильин — Новосибирск: Изд-во ИВМиМГ СО РАН, 2007. — 371 с.
5. Ильин, В.П. Параллельные методы и технологии декомпозиции областей / В.П. Ильин // Вестник ЮУрГУ. Серия «Вычислительная математика и информатика». — 2012. — No. 46(305). — С. 31-44.
6. Ильин, В.П. Параллельные методы декомпозиции в пространствах следов / В.П. Ильин, Д.В. Кныш // Вычислительные методы и программирование. — 2011. — Т. 12, No. 1.
— С. 100-109.
7. Saad, Y. Iterative Methods for Sparse Linear Systems, Second Edition / Y. Saad — SIAM, 2003. — 528 p.
8. Intel Math Kernel Library. Reference Manual. URL: http://software.intel.com/ sites/products/documentation/doclib/mkl_sa/11/mklman/index.htm (дата обращения:
12.02.2013).
102 Вестник ЮУрГУ. Серия «Вычислительная математика и информатика»
9. Ильин, В.П. Методы бисопряженных направлений в подпространствах Крылова / В.П. Ильин // СибЖИМ. — 2008. — Т. 11, No. 4(36). — С. 47-60.
10. Bell, N. Cusp: Generic Parallel Algorithms for Sparse Matrix and Graph Computations / N. Bell, M. Garland. URL: http://cusp-library.googlecode.com (дата обращения:
15.10.2012).
11. PETSc: Home Page. URL: http://www.mcs.anl.gov/petsc/ (дата обращения:
12.02.2013).
12. Hypre. URL: http://acts.nersc.gov/hypre/ (дата обращения: 12.02.2013).
13. Yousef Saad - Software. URL: http://www-users.cs.umn.edu/~saad/software/ (дата обращения: 12.02.2013).
14. Кластер НКС-30Т. URL: http://www2.sscc.ru/HKC-30T/HKC-30T.htm (дата обращения:
12.02.2013).
15. Nabben, R. A comparison of abstract versions of deflation, balancing and additive coarse grid correction preconditioners / R. Nabben, C. Vuik // Numerical Linear Algebra with Applications. — 2008. — Vol. 15, No. 4. — P. 355-372.
Дмитрий Сергеевич Бутюгин, младший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН, аспирант, Новосибирский государственный университет (Новосибирск, Российская Федерация), [email protected].
Яна Леонидовна Гурьева, к.ф.-м.н., старший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), [email protected].
Валерий Павлович Ильин, д.ф.-м.н., профессор, главный научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), [email protected].
Данил Валерьевич Перевозкин, младший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), [email protected].
Артем Владимирович Петухов, младший научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), [email protected].
Игорь Николавевич Скопин, к.ф.-м.н., научный сотрудник, Институт вычислительной математики и математической геофизики СО РАН (Новосибирск, Российская Федерация), [email protected].
PARALLEL ALGEBRAIC SOLVERS LIBRARY KRYLOV
D.S. Butyugin, Institute of Computational Mathematics and Mathematical
Geophysics SB RAS (Novosibirsk, Russian Federation),
Y.L. Guryeva, 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),
A.V. Petukhov, Institute of Computational Mathematics and Mathematical
Geophysics SB RAS (Novosibirsk, Russian Federation),
I.N. Skopin, Institute of Computational Mathematics and Mathematical Geophysics
SB RAS (Novosibirsk, Russian Federation)
Article describes functional capabilities and software implementation peculiarities of parallel algorithms library Krylov, which is oriented on the solution of large systems of linear algebraic equations with sparse symmetric and unsymmetric matrices (positive definite and semi-definite) obtained from discrete approximations of multidimensional boundary value problems for partial differential equations on unstructured meshes. The library includes two-level iterative methods in Krylov subspaces; preconditioning of the latter is based on the balanced decomposition of the computational domain with variable sizes of subdomain overlapping and different boundary conditions on interfacing boundaries. Program implementations use typical compressed sparse matrix data formats. Results of numerical experiments are presented which demonstrate the efficiency of parallelization for typical ill-conditioned problems.
Keywords: preconditioned iterative algorithms, Krylov subspaces, domain decomposition methods, sparse algebraic systems, numerical experiments.
References
1. Butyugin D.S., Il’in V.P., Itskovich Y.A., et al. Krylov: biblioteka algoritmov i programm dlya resheniya SLAU [Krylov: library of algorithms and programs for SLAEs solution]. Sovremennye problemy matematicheskogo modelirovaniya. Matematicheskoe modelirovaniye, chislennye metody i compleksy programm. Sbornik trudov Vserossiyskikh nauchnykh molodezhnykh shkol [Modern problems of mathematical simulation. Mathematical simulation, numerical methods and program complexes. All-Russian young scientists schools proceedings]. Rostov-na-Donu, Publishing of South Federal University, 2009. P. 110-128.
2. Butyugin D.S. Metody parallelnogo resheniya SLAU na systemakh s raspredelennoy pamatyu v biblioteke Krylov [Methods of parallel SLAEs solution on the systems with distributed memory in Krylov library]. Vestnik YUURGU. Seriya “Vychislitelnaya matematika i informatika” [Bulletin of South Ural State University. Seriers: Computational Mathematics and Software Engineering], 2012. No. 47(306). P. 5-19.
3. Il’in V.P. Metody konechnykh raznostey i konechnykh objemov dlya ellipticheskikh uravnenij [Methods of finite differences and finite volumes for elliptic equations] [Methods and technologies of finite elements]. Novosibirsk, ICM&MG SBRAS Publishing, 2001.
4. Il’in V.P. Metody i tekhnologii konechnykh elementov [Methods and technologies of finite elements]. Novosibirsk, ICM&MG SBRAS Publishing, 2007.
5. Il’in V.P. Parallelnye metody i tekhnologii decompozitsii oblastey [Parallel methods and technologies of domain decomposition]. Vestnik YUURGU. Seriya “Vychislitelnaya matematika i informatika” [Bulletin of South Ural State University. Seriers: Computational Mathematics and Software Engineering], 2012. No. 46(305). P. 31-44.
6. Il’in V.P., Knysh D.V. Parallelnye metody decompozitsii v prostranstvakh sledov [Parallel decomposition methods in trace spaces]. Vychislitelnye metody i programmirovanie [Computational methods and programming], 2011. Vol. 12, No. 1. P. 100-109.
7. Saad Y. Iterative Methods for Sparse Linear Systems, Second Edition. SIAM, 2003.
8. Intel Math Kernel Library. Reference Manual. URL: http://software.intel.com/
sites/products/documentation/doclib/mkl_sa/11/mklman/index.htm (accessed:
12.02.2013).
9. Il’in V.P. Metody bisopryazhennykh napravleniy v podprostranstvakh Krylova [Bi-conjugate directions methods in Krylov subspaces]. SibZHIM [Syberian Journal of Industrial Mathematics], 2008. Vol. 11, No. 4(36). P. 47-60.
10. Bell, N. Cusp: Generic Parallel Algorithms for Sparse Matrix and Graph Computations j N. Bell, M. Garland. URL: http://cusp-library.googlecode.com (accessed: 15.10.2012).
11. PETSc: Home Page. URL: http://www.mcs.anl.gov/petsc/ (accessed: 12.02.2013).
12. Hypre. URL: http://acts.nersc.gov/hypre/ (accessed: 12.02.2013).
13. Yousef Saad - Software. URL: http://www-users.cs.umn.edu/~saad/software/ (accessed:
12.02.2013).
14. Klaster HKC-30T [Cluster HKC-30T].
URL: http://www2.sscc.ru/HKC-30T/HKC-30T.htm (accessed: 12.02.2013).
15. Nabben R., Vuik C. A comparison of abstract versions of deflation, balancing and additive coarse grid correction preconditioners jj Numerical Linear Algebra with Applicati 2008. Vol.
15, No 4. P. 355-372.
Поступила в редакцию 14 июня 2013 г.
