CC BY
114
УДК 519.6
0 ВОПРОСАХ РАСПАРАЛЛЕЛИВАНИЯ КРЫЛОВСКИХ ИТЕРАЦИОННЫХ МЕТОДОВ1
В.П. Ильин
В работе рассматриваются математические вопросы многообразных вычислительных технологий методов распараллеливания итерационных процессов крыловского типа для решения больших разреженных симметричных и несимметричных СЛАУ, возникающих при сеточных аппроксимациях многомерных краевых задач для систем дифференциальных уравнений. Характерным примером являются конечно-элементные приближения в газогидродинамических приложениях, где в каждом узле определены пять неизвестных функций, в силу чего СЛАУ имеет мелкоблочную структуру. Основой применяемых алгоритмов является гибкий метод обобщенных минимальных невязок FGMR.ES с динамическими пре-добуславливателями аддитивного типа, представляющий собой верхний уровень двухступенчатого итерационного алгоритма Шварца.
Для повышения производительности алгебраических решателей автором предлагается применение различных подходов: декомпозиции расчетной области с различными топологиями, типами краевых условий на смежных границах и размерами пересечений подобластей, методов грубосеточной коррекции и агрегации, дефляции и неполной факторизации матриц. Описываются унифицированные формулировки используемых алгоритмов, а также вопросы их вычислительной эффективности и масштабируемого распараллеливания на суперкомпьютерах гетерогенной архитектуры. Приводятся примеры технологических требований к особенностям программных реализаций библиотек параллельных алгоритмов для решения систем линейных алгебраических уравнений.
Ключевые слова: итерационные методы, подпространства Крылова, предобусловленные матрицы, декомпозиция областей, параллельные алгоритмы, программные и вычислительные технологии.
Введение
В последние годы сформировалось бурно развивающееся направление вычислительной алгебры, связанное с решением больших разреженных систем линейных алгебраических уравнений (СЛАУ), возникающих из сеточных аппроксимаций (конечно-элементных или конечно-объемных, см. [1] и цитируемые там работы) многомерных краевых задач для систем дифференциальных уравнений в частных производных. Естественно, это связано с актуальными проблемами математического моделирования для сложных междисциплинарных приложений, включающих гидро-газодинамические процессы, динамику напряженно-деформированных состояний и др., на современных архитектурах суперкомпьютеров пе-тафлопных масштабов. Само понятие большой задачи быстро эволюционирует, и сейчас необходимо говорить о решениях СЛАУ с порядками 1010 на многопроцессорных вычислительных системах (МВС) с числом ядер, или потоков, в десятки и сотни тысяч.
Повышение производительности вычислений, или быстродействия, в рассматриваемой прикладной сфере обуславливается двумя основными факторами. Первый — это конструирование итерационных процессов с высокой скоростью сходимости (прямые алгоритмы применимы только для относительно небольших размерностей задач), а второй — технологическое обеспечение масшатабируемости распараллеливания, что означает сохранение
1 Статья рекомендована к публикации программным комитетом Международной научной конференции «Параллельные вычислительные технологии 2013».
эффективности с ростом порядков СЛАУ и/или количества вычислительных устройств. Архитектура последних развивается главным образом количественно, а не качественно, и типовая структура гетерогенной МВС — это кластер с одинаковыми или почти одинаковыми вычислительными узлами, каждый из которых содержит несколько центральных процессорных многоядерных элементов (CPU) с общей памятью, а также несколько специализированных ускорителей, среди которых наибольшее распространение получили «графические процессорные элементы общего назначения» (GPGPU, или просто GPU) с очень большим числом ядер и сложной организацией внутренней иерархической оперативной памяти. Реализация параллельных вычислений на такого типа МВС осуществляется средствами систем MPI, OpenMP и CUDA, причем оптимизация программного кода фактически производится разработчиком в основном вручную, поскольку соответствующие средства автоматизации недостаточно функциональны.
Главными инструментами для решения возникающих здесь алгоритмических проблем являются методы декомпозиции областей [2-8] и итерационные процессы в подпространствах Крылова [5, 1]. Хотя эти направления уже стали классическими в вычислительной математике, они сейчас переживают период бурного развития, и буквально каждый год появляются новые интересные идеи. Рост эффективности алгебраических решателей достигается с помощью применения многообразных подходов: оптимизация размеров и топологий пересечений соседних подобластей, а также типов краевых условий на смежных границах, методы грубосеточной коррекции и агрегации, дефляции и неполной факторизации.
В разделе 1 мы даем постановку и особенности рассматриваемых задач, а также сравнительный анализ методов их решения. Раздел 2 посвящен как общим, так и специальным технологическим аспектам распараллеливания алгоритмов, а в последнем разделе описываются примеры технологических требований к программным реализациям алгебраических решателей и, в частности, к библиотеке KRYLOV, предварительные сообщения о которой опубликованы в [6, 8]. В заключении отмечается комплексность взаимосвязанных теоретических, алгоритмических и технологических вопросов, определяющих уровень высокопроизводительных вычислений математического и программного обеспечения для решения актуальных алгебраических систем.
1. Математические аспекты двухуровневых методов декомпозиции
Рассматривается следующая задача: пусть вещественная матрица СЛАУ
разбита на блочные строки А = {Ар є 'Я,Мр’М, р = 1, ...,Р}, N1 + ... + N = N с приблизительно равным числом строк N ^ 1 в каждой. Соответствующим образом также разбиваются векторы искомого решения и правой части и = {ир}, / = {/р}, так что система уравнений (1) может быть записана в форме совокупности Р подсистем
Au = f, A = {aij} є Rn,n,
(1)
P
(2)
q = 1
где Ар,д є — прямоугольные матричные блоки, получаемые при разбиении каждой
блочной строки (и всей матрицы А) на блочные столбцы.
Пусть матрица (1) задана в сжатом разреженном СБИ-формате [3] с указанием числа строк N в каждой из Р подсистем (2), в соответствии с разбиением матрицы А на блочные строки Ар. Естественно предполагается, что строки исходной матрицы пронумерованы
р-1 р
подряд от 1 до N, так что номера строк каждого блока Ар меняются от 1 + ^ N до ^ N
г=1 г=1
(эти номера назовем глобальными).
Предполагается, что СЛАУ (1) представляет собой систему сеточных уравнений, так что каждая компонента векторов и, / соответствует узлу сетки, общее число которых в
р
расчетной сеточной области = У 0р равно N. При этом блочное разбиение матриц и
р=1
векторов соответствует разбиению (декомпозиции) на Р сеточных непересекающихся подобластей 0р, в каждой из которых находится Np узлов.
Описанная декомпозиция не использует узлов — разделителей, т.е. условные границы подобластей не проходят через узлы сетки. Формальности ради можно считать, что внешность расчетной области есть неограниченная подобласть 0о с числом узлов N0 = 0. Сеточное уравнение для г-го узла сетки может быть записано в виде
ai,iui 1 ^ ' ai,juj — г € 0 , (3)
где через Шi обозначается сеточный шаблон г-го узла, т.е. совокупность номеров всех узлов, участвующих в г-м уравнении.
Будем считать, что сеточные узлы и соответствующие переменные пронумерованы следующим образом: сначала идут подряд все узлы 1-й подобласти О1 (неважно в каком внутреннем порядке), затем — узлы второй подобласти и т.д. Отметим, что взаимнооднозначного соответствия алгебраической и геометрической (сеточной) интерпретации структуры СЛАУ может и не существовать. Типичный пример — одному узлу сетки может соответствовать т > 1 переменных в алгебраической системе. Если для каждого узла такие кратные переменные нумеруются подряд, то структура СЛАУ приобретает мелкоблочный (т ^ N) вид (в (3) щ и /i означают не числа, а подвекторы порядка т соответственно, € ^т’т), и на таких случаях мы остановимся в последующем особо.
1.1. О построении параллельных предобуславливателей
Для заданного блочного разбиения СЛАУ, или декомпозиции сеточной расчетной области без пересечения подобластей, можно построить хорошо распараллеливаемый аддитивный метод Шварца, представимый итерационным процессом «по подобластям»:
р
Ар,рир = /р — Ар>4 1 = Ур 1- (4)
ч=1
Здесь п будем называть номером внешней итерации, которая в действительности реализуется более утонченным образом, чему будет посвящен специальный раздел.
В целом алгоритм решения крупноблочной СЛАУ является двухуровневым и нижний, или внутренний, уровень заключается в решении (прямым или итерационным методом) для каждого п — независимо и параллельно — подсистем вида (4) с матрицами Ар,р.
При вычисленных правых частях урЩ-1 нахождение всех требует для каждой подобласти решения подсистем с порядками ^, которые на каждой п-й итерации могут выполнять-
ся независимо и одновременно на разных процессорах. В данном случае перевычисление урП-1 фактически означает использование новых значений и^-1 для соседних подобластей в качестве граничного условия 1-го рода (Дирихле) для околограничных внутренних узлов из 0р.
Известно, что скорость сходимости итераций метода Шварца (4) возрастет, если декомпозицию расчетной области сделать с пересечением подобластей.
В силу этого мы дадим определение расширенной сеточной подобласти 0р Э 0р, имеющей пересечения с соседними подобластями, величину которых мы будем определять в терминах количества околограничных сеточных слоев, или фронтов. Обозначим через Гр € 0р множество внутренних околограничных узлов из 0р, т.е. таких узлов Р^ € 0р, у которых один из соседей не лежит в 0р(Р?- € 0р, ] € ] — г)). Обозначим далее через Г1 мно-
жество узлов, соседних с узлами из Гр, но не принадлежащих 0р, через Гр - множество узлов, соседних с узлами из но не принадлежащих объединению Гр и 0р, через Гр -множество узлов, соседних с Г^, но не принадлежащих Г У Гр У 0р, и т.д. Соответственно эти множества назовем первым внешним слоем (фронтом) узлов, вторым слоем, третьим и т.д. Получаемое объединение узлов
0р = при Гр-и г^
будем называть расширенной р-й сеточной подобластью, а целую величину А определяем как величину расширения, или пересечения (в терминах количества сеточных слоев). Случай А — 0 фактически означает декомпозицию области 0^ на подобласти без пересечений (00 = 0р).
На формальном алгебраическом языке каждой расширенной подобласти можно сопоставить подсистему уравнений
р
(Ар,р + в^р)ир = ,/р ^ ' Ар,9+ в^рир (5)
4 = 1 4=Р
где ^р — диагональная матрица, определяемая соотношением
р
А>е = ^2 Ар>д е, е = (1,..., 1)Т €ПЯр.
4 = 1 4 = Р
Здесь размерность векторов —П и /р равна числу узлов N в расширенной подобласти Ар,р, ^р € , а введенный параметр в формально позволяет рассматривать зависящие
от ир и — правые части (5) как аппроксимации краевых условий на смежных границах
подобластей: в = 0 соответствует условию 1-го рода (Дирихле), в = 1 — условию 2-го рода
(Неймана), а в € (0,1) — условию 3-го рода (Робена), [10].
Переходя теперь к полному вектору и и вводя блочно-диагональные матрицы
Вр = Ыосй-^гау{Ар’р + в^р}, (6)
из (5) получаем итерационный процесс вида
ип = ип-1 + В-1(/ - Аип-1) = ип-1 + В-1гп-1 =
- - (7)
= Тип-1 + / , Т = I - В-1А, / = В-1/,
где B — предобуславливающая матрица, определяемая следующим образом.
Введем матрицу Wp = (wi,..., wnp)T £ RNp’N “сужения” вектора u £ RN в пространство RWp соответствующей подобласти, для чего компоненты каждой вектор-строки wk = {wfc’i} полагаются равными единице, если i £ Qp, и нулю в противном случае. При этом Wp^ является матрицей оператора продолжения пространства Qp в Q, и в итоге предобуславливатель аддитивного метода Шварца (additive Swartz) принимает вид
P
B-1 = B-S = Е Bp-1, Bp-1 = WpTBp-1Wp. (8)
p=i
На алгебраическом языке формулы (7) при B = Bas определяют стационарный блочный метод Якоби (BJ), на каждом шаге которого надо обращать матрицу Bp, что фактически означает одновременное (параллельное) решение расширенных СЛАУ в подобластях. Если эти процедуры осуществляются с помощью “внутренних” итерационных процессов, то это неизбежно приводит к тому, что в (7) на каждом n-м шаге реально используется переменное (динамическое) предобуславливание с матрицами Bn.
Кардинальный подход к ускорению итераций заключается в переходе от формул (7) для BJ к алгоритмам в предобусловленных подпространствах Крылова Kn(v0, v1,..., vn-1), где введены обозначения v0 = B0"1ro,vfc = B-1Avk-1. Один из возможных путей здесь заключается в применении гибкого (flexible) метода обобщенных минимальных невязок FGMRES [5], который на каждом n-м шаге минимизирует норму ||rn||, но требует для своей реализации оперативную память объемом 2nN. Использование FGMRES с рестартами через каждые m итераций при больших n позволяет значительно сокращать объем требуемой памяти до 2mN, но существенно замедляет скорость сходимости.
Главный, и практически единственный, путь повышения эффективности алгоритмов — это развитие крыловских методов или за счет совершенствования предобуславливателей, или путем улучшения итерационных подпространств (добавление новых базисных векторов, в частности, что тоже может интерпретироваться как введение дополнительного предобу-славливания).
Одним из резервов ускорения является грубосеточная коррекция (coarse grid correction, CGC), или агрегирование, а также идейно примыкающие сюда алгебраические многосеточные методы (algebraic multi-grid, AMG). Суть здесь заключается в том, что на каждом шаге стандартного BJ эволюция последовательных приближений от каждой области передается только ее непосредственным соседям. Идея исправления данной ситуации — формирование и решение дополнительных относительно небольших СЛАУ, которые соединяли бы все подобласти и передавали бы, пусть грубо, глобальные итерационные возмущения.
Реализация данного подхода заключается в конструировании дополнительного предо-буславливателя и внешне аналогична (8). Мы формируем новый оператор сужения с матрицей Wc £ RNc,N, Nc ^ N, в которой каждая строка содержит только по несколько ненулевых элементов (обычно равных единице), соответствующих одной из подобластей. Транспонированная матрица W^ будет представлять тогда оператор продолжения, грубосеточный предобуславливатель принимает вид
B-1 = WcTA-1WC, Ac = WcAWcT £ RNc’Nc, (9)
и полный предобуславливатель определяется аддитивным образом:
в-1 = в-£ + Bc-1.
(10)
Отметим, что реализация СОС на каждой п-й итерации требует дополнительного решения СЛАУ с матрицей невысокого порядка Ас, что можно сделать с помощью прямого решателя, например, РАКБШО из библиотеки МКЬ ШТЕЬ [9].
Развитие идей грубосеточной коррекции активно продолжается в различных направлениях, см. [11, 12] и цитируемые там работы. Многосеточные методы используют технологии интерполяции на каждом из последовательных этапов дискретизации. Адаптивные алгоритмы сглаженного агрегирования основаны на взвешенных усреднениях различных приближений. Методы типа РЕТ1 применяют декомпозицию области с явным выделением разделителей сеточных подмножеств (макрограни, макроребра, вершины) и построение дополнений Шура на иерархическом принципе.
Опишем еще один подход к ускорению крыловских процессов — метод дефляции, который имеет широкое распространение в разных версиях. Мы его представим в применении к выбору начального приближения, если СЛАУ с одинаковой матрицей решается многократно с разными правыми частями. Именно такая ситуация возникает в двухуровневых методах декомпозиции областей. Пусть в результате предыдущего решения СЛАУ с помощью какого-то крыловского метода вычислены А-ортогональные векторы (^1,..., -шт) = ,
составляющие базис пространства, которое будем называть дефляционным. Для решения новой системы выбираем начальное приближение и0 таким, чтобы соответствующая начальная невязка г0 и начальный направляющий вектор р0 удовлетворяли условиям ортогональности
Допустим, что и 1 есть «предварительный» начальный вектор, а г 1 = / — Аи 1 — соответствующая невязка. Тогда условия (11) будут выполняться, если положить
В рамках одной статьи невозможно дать содержательный обзор современных тенденций в развитии предобусловленных итерационных процессов крыловского типа. Можно только констатировать, что оригинальные подходы появляются практически непрерывно, и данный момент следует рассматривать как важный технологический фактор при создании математического и программного обеспечения, ориентированного на длительный жизненный цикл.
В заключение данного раздела отметим, что если в матрице исходной СЛАУ (1) выделить главную блочную диагональ (А = В — С, В = {Ар,р}), то при условии ее разложения на блочно-треугольные множители систему можно переписать в следующем виде:
"Тг0 = 0, "ТАр0 = 0.
(11)
и0 = и-1 + "А-1"Тг-1, г0 = / — Аи0, р0 = [I — "А-1(А"Т)]г0, АЙ = "ТА".
(12)
Ви = Си + /, В = ^= (£д^), и = Ьд/Сид1^ и + Ьд/.
(13)
Отсюда можно ввести формально предобусловленную СЛАУ
Аи = (I — Т)и = / , и = Ц/и,
Т = ьд1сид1(= Т*), / = ьд1/,
(14)
для решения которой естественно применять блочный (двусторонне предобусловленный) метод Якоби
un = Tun-1 + f, un = Udun. (15)
Заметим, что если матрица A симметрична, то этим же свойством обладают и матрицы T, A. Отметим, что в формулах (13) - (15) можно положить
Ld = D, Ud = I, или Ld = I, Ud = D,
что будет соответствовать одностороннему предобуславливанию (левому или правому соответственно).
1.2. Некоторые особенности методов в подпространствах Крылова
Как видно из предыдущего раздела, итерационные методы (мы подразумеваем «в подпространствах Крылова», как оптимальные) можно применять для решения как «глобальной» СЛАУ, т.е. для ускорения блочного метода Якоби, так и для «локальных» систем в подобластях. В дальнейшем крыловские алгоритмы мы рассмотрим в единообразной форме, независимо от того, применяются они на внешних или внутренних итерациях, основой которой является А^-ортогонализация различных векторов, s = 0 или s = 1, см. [1].
Если предобусловленная СЛАУ Au = f симметрична, то наиболее экономичными являются методы сопряженных градиентов или сопряженных невязок (CG или CR, для s = 0 или s = 1 соответственно), описываемые следующими формулами:
г0 = f — Au0 = u1 — u0, u1 = Tu0 + f, p0 = r0,
un+1 = un + a4s)pn, a4s) = pis)/^is), pis) = (Asrn, rn), d^)
^ = (Apn, A(s)pn), rn+1 = rn — ais)Apn, ( )
pn+1 = rn+1 + ^V, eis) = piS+1/pis).
Различные методы решения несимметричных систем базируются на As-
ортогонализация Арнольди с предобуславливанием. Поскольку в общем случае предобу-славливающие матрицы Bn могут отличаться на разных итерациях, то мы рассматриваем «гибкую» (flexible) ортогонализацию в следующем виде:
un = u0 + y1V1 + ... + y„vn, (vn, Asvk) = dis)4’„,
dis) = (vn, Asvn),
vn+1 = AB-1vn — ЕП=1 , v1 = r0 = f — Au0, n = 1,2,...
, (s) = (Avn’A°vfc) k = 1 n + 1 V , = (v1 vn+1) (17)
= (A=vfc,vk) , k = 1,...,n + 1, Vn+1 = (v ,...,v )
—{hfc,ra}—
e*
є Rra+1>ra, Яга є
Получаемые предобусловленные векторы Bfc 1v^ образуют базис подпространства Кры-
лова
Кга+1(г0, А) = 8рап(В- V,..., Вга+ ^га+1} = 8рап(В- 1г0, АВ2 1г0,..., А^В-^г0}, (18)
и на их основе формируются алгоритмы полной А^-ортогонализации или методы обобщенных А^-минимальных невязок (РОМ, или А-РОМ, ОМИЕЯ, или А-ОМИЕЯ), см. [1, 5].
54 Вестник ЮУрГУ. Серия «Вычислительная математика и информатика»
Внешние итерации оканчиваются по выполнению условия
(rn,rn) < 4f/),
где eex ^ 1 характеризует точность решения. Оптимизация критериев останова внутренних итераций — это непростой вопрос, и здесь параметры точности еП в принципе должны зависеть от номера n внешней итерации.
2. Технологические вопросы реализации параллельных итерационных алгоритмов
Рассматриваемые в данном разделе вопросы — это конкретизация общей проблемы «отображение алгоритмов на архитектуру МВС». В качестве типовой модели вычислительной системы может служить НКС-30Т — гетерогенный кластер ИВМиМГ СО РАН [13]. Такая МВС формально представляет набор вычислительных узлов, которые с помощью коммуникационной сети обеспечиваются связями «каждый с каждым», управляемыми с помощью программных средств системы передачи сообщений MPI [14]. Каждый узел имеет свою многоуровневую память (большую оперативную и поменьше — сверхбыстрый кэш, тоже имеющий внутренние уровни), общую для расположенных в нем многоядерных центральных процессорных устройств (CPU, в НКС-30Т их два, каждый с четырьмя или шестью ядрами), а также графические процессорные устройства с очень большим числом ядер (GPU, в НКС-30Т их три, с 512 ядрами в каждом). На каждом CPU (и даже на их отдельных ядрах) можно формировать вычислительный MPI-процесс, внутри которого распараллеливание реализуется средствами системы OpenMP [15] над общей памятью (более конкретно — организуется несколько вычислительных потоков). На одном GPU допускается запуск только одного MPI-процесса, причем внутри него программирование осуществляется на языке CUDA [16] с достаточно сложными средствами управления внутренней иерархической памятью. Особенностью (и недостатком) GPU является относительно медленные коммуникации с памятью CPU. В целом средства MPI обеспечивают организацию синхронных или асинхронных вычислений, в том числе при совмещении их во времени с передачей данных от процессора к процессору.
Как видно из приведенного описания, построение даже грубой математической модели вычислительного процесса на такой архитектуре, с целью его оптимизации и оценок коммуникационных потерь, не представляется возможным. Поэтому соответствующие исследования являются сугубо экспериментальными, а основной инструмент в данном случае — это метод проб и ошибок.
Справедливости ради следует сказать, что математик-программист работает не с пустыми руками, а при наличии достаточного богатого вычислительного инструментария (главным образом библиотеки BLAS и SPARSE BLAS [8], содержащие основные алгебраические операции с векторами и матрицами, в том числе разреженными), созданного профессионалами с помощью экономичных языковых средств низшего уровня. В частности, можно упомянуть библиотеку CUSP [17] для решения на GPU разреженных СЛАУ с использованием средств BLAS.
По поводу вопросов реализации параллельных итерационных алгоритмов мы отметим две проблемы, относящиеся к области вычислительных технологий. Первая из них касается формирования вспомогательных СЛАУ для подобластей, которые необходимо решать
на нижнем уровне двухуровневых итерационных методов декомпозиции областей. В силу существующего для большинства случаев изоморфизма сеточных шаблонов и портретов матричных строк анализ декомпозиции может проводиться одинаковым образом как в терминах сеточных графов, так и на языке матричных графов, в силу чего употребляются названия «алгебраическая декомпозиция» и «геометрическая декомпозиция». Более целесообразным и естественным выглядит реализация декомпозиции на этапе построения сетки, когда можно оперировать информацией о геометрических объектах расчетной области, топологическими связями и расстояниями, и т.д. Однако на практике зачастую при использовании библиотеки алгебраических решателей декомпозиция не задается, и в таких случаях фактически требуется построить по каким-то критериям блочную структуру СЛАУ, пользуясь только матричным CSR-форматом. Эта задача имеет достаточно высокую информационную и логическую сложность, особенно если ее требуется решить в параллельном, т.е. распределенном по различным процессорам, режиме.
Формирование пересекающихся сеточных подобластей и соответствующих расширенных СЛАУ осуществляется в два этапа. На первом определяются сбалансированные (с примерно равным числом узлов) подобласти без перехлеста, для чего могут использоваться инструментарии типа популярной библиотеки METIS [18] (в том числе существующей в параллельном варианте), осуществляющей операции над графами. При этом подобласти каким-либо образом упорядочиваются (это соответствует разбиению матрицы на блочные строки), и в соответствии с этим нумеруются узлы сетки (или матричные строки): сначала идут все узлы (строки) из первой подобласти (блочной строки), затем — из второй, и т.д. На втором этапе производится постепенное расширение подобластей по слоям соседних узлов, или фронтов: сначала к внутренним узлам подобласти присоединяются ближайшие внешние узлы, являющиеся непосредственными соседями к внутренним и которые образуют 1-й слой расширения, затем к ним аналогично присоединяются узлы 2-го слоя, и т.д. до заданного априори количества фронтов расширения.
В матричной технологии это осуществляется только на основе глобального CSR-формата, а результатом являются локальные матричные CSR-форматы для каждой из расширенной подобластей (при этом, естественно, требуется перенумерация матричных строк и столбцов соответствующих ненулевых элементов из глобальной упорядоченности в локальную). Если матрицы расширенных СЛАУ в подобластях формируются только один раз, то их правые части необходимо пересчитывать в околограничных узлах на каждой внешней итерации. Для этого требуется определять совокупности множеств узлов-доноров и узлов-акцепторов, которые передают информацию от своей подобласти к соседним и, наоборот, принимают данные.
Второй требующий решения технологический вопрос — как и в каких пространствах осуществлять внешний итерационный процесс? Наиболее прямой и естественный путь заключается в реализации крыловского метода на основе формулы (7), где un и f полные векторы с размерностью, равной порядку исходной СЛАУ (1), т.е. O(h-3) в трехмерном случае, где h есть характерный шаг сетки. В этом случае метод FGMRES исполняется в кластерном варианте с помощью P MPI-процессов, что требует, в частности, дополнительных обменов при вычислении скалярных произведений векторов.
Однако существует и альтернативный подход, состоящий в проведении внешних итераций в пространстве следов, т.е. для векторов, определенных только на смежных границах пересекающихся подобластей [10]. В этом случае исходная СЛАУ формально редуциру-
ется путем исключения неизвестных, соответствующих внутренним узлам подобластей, и итоговая размерность системы понижается на порядок, т.е. до O(h-2). При этом реализацию внешнего FGMRES можно осуществлять только на одном процессоре, используя распараллеливание ограниченными средствами OpenMP, но полностью избегая при этом коммуникационных потерь. Однако такая вычислительная технология имеет существенный недостаток — неизбежный дефицит оперативной памяти одного процессора при решении больших СЛАУ с масштабируемым параллелизмом, требующим формирования очень большого количества подобластей.
Что касается не вычислительных, а программных и информационных технологий ускорения параллельных процессов для рассматриваемых классов задач и алгоритмов, то здесь также имеются значительные резервы за счет оптимизации кода, организации развертки циклов, выбора режимов компиляции, профессионального использования командных возможностей систем MPI и OpenMP, и т.д. Конечный результат здесь в существенной степени зависит от искусства математика-программиста в синхронизации массовых расчетных и обменных операций, имея целью добиться максимум возможного в условиях ограниченных рамок маневрирования вычислительными потоками на фиксированных функциональных характеристиках аппаратных устройств МВС. Но так или иначе, человеческий фактор в существующих условиях еще играет далеко не последнюю роль.
3. Примеры технических требований к библиотекам алгебраических решателей
Программного обеспечения по вычислительной алгебре в мире существует достаточно много и давно, как самостоятельного, так и в составе пакетов прикладных программ, как коммерческого, так и общедоступного. Однако библиотек по итерационным параллельным алгоритмам решения больших разреженных СЛАУ имеется не так уж много. Помимо уже упоминавшейся библиотеки CUSP для GPU, можно назвать такие разработки, как Hypre [19], PETSc [20] и Saad Software [21]. Высокопроизводительные продукты для решения сверхбольших СЛАУ на суперкомпьютерах петафлопного масштаба еще ждут своего часа, так что мы рассмотрим некоторые технические требования к библиотеке алгебраических решателей, руководствуясь в первую очередь практическими требованиями, возникающими из актуальных проблем математического моделирования. Более того, мы будем исходить из целевой постановки о создании программного обеспечения, ориентированного на широкого круга пользователей и рассчитанного на длительный жизненный цикл с регулярным развитием функционального наполнения, а также адаптирующегося к эволюции компьютерных архитектур.
а. Состав решаемых задач и алгоритмов. Обладая определенным багажом знаний и технологий, или ноу-хау, заманчиво на единых принципах охватить СЛАУ различных типов: вещественных и комплексных, эрмитовых и неэрмитовых, положительно определенных и знаконеопределенных и т.д. Такую систематизацию можно продолжить до выделения специальных систем, для которых применимы сверхбыстрые алгоритмы. Основа современных итерационных процессов — это предобусловленные методы в крыловских подпространствах. Однако в этих обоих направлениях имеется огромное разнообразие вариантов, и выше мы только коротко остановились на двухуровневых методах декомпозиции областей.
б. Принципы организации интерфейса. Традиционный вопрос в данном случае такой — делать ли алгебраический решатель в форме «черного ящика» или, наоборот, в качестве полностью открытого пользователю и настраиваемого на конкретные задачи инструмента? Возможны, конечно, и различные компромиссные «серые» варианты. Целесообразные решения, естественно, зависят от того, на какого типа пользователя рассчитан программный продукт. Наиболее реальной является ситуация, когда решатель нужен не в автономном варианте, а встраиваемый в уже существующий ППП. Традиционная форма универсальной программной реализации итерационного метода крыловского типа такова: широкое использование BLAS-овских функций для векторных операций, а также открытость для внешних процедур умножения вектора на матрицу исходной СЛАУ и на предобуславливатель. Такое абстрагирование изначально закладывается в концепцию объектно-ориентированного программирования, например, в языке С++. Однако хорошо известно, что чрезмерный универсализм — это враг эффективности и экономичности.
Отсюда возникает немаловажный вопрос — как встраивать в широко-форматную библиотеку специальные сверхбыстрые алгоритмы для частного вида СЛАУ, которые не укладываются в общую вышеприведенную схему? Например, расчетная область или даже одна из подобластей могут допускать разделение переменных и, как следствие, применение быстрого преобразования Фурье (БПФ), которого очень жалко было бы запрещать к использованию.
в. Проблема переиспользования программ. Разработанное компьютерным сообществом прикладное программное обеспечение представляет огромный интеллектуальный потенциал, и умение его использовать дает заведомое преимущество разработчикам программного продукта. В рассматриваемых нами обстоятельствах это означает, в первую очередь, наличие первоклассных программ, реализующих прямые методы, которые можно и нужно использовать для решения СЛАУ в подобластях.
г. Внутренняя структура и организация библиотеки. Один из возможных способов построения библиотеки алгебраических решателей — это создание каталогизированного и систематизированного хранилища большого количества алгоритмов, из которого специальными конфигурационными средствами может формироваться конкретная версия продукта для определенных условий эксплуатации. Более примитивный подход — конкретная программа просто вынимается из хранилища (переписывается в файл) и отторгается. Третий, и наиболее продвинутый, способ существования прикладной библиотеки состоит в создании баз данных для типовых СЛАУ, для результатов их решения различными алгоритмами, а также в организации тренинга и обучения потенциальных пользователей.
Список технологических и примыкающих организационных вопросов по созданию библиотеки параллельных алгоритмов решения больших разреженных СЛАУ можно было бы значительно продолжить (многоязычность и многоверсионность, платформонезависимость, внутренние средства развития и адаптируемости, сопровождение, документируемость, ли-цензионность и т.п.). Частично эти вопросы были решены разработчиками библиотеки KRYLOV, которая находится на стадии опытной эксплуатации в ИВМиМГ СО РАН в авторском сопровождении.
Режим исполнения в библиотеке KRYLOV можно назвать полуавтоматическим, или в стиле «серого ящика». Внешний итерационный процесс по входному заданию пользователя
может выполняться или на корневом процессоре в пространстве следов, или в распределенном варианте. Для решения внутренних СЛАУ в подобластях допустимо использование как авторских реализаций различных крыловских процессов, так и прямой алгоритм PARDISO, а также ряд решателей из библиотеки CUSP NVIDIA. В распоряжении пользователя - различные методы декомпозиции с параметризованными размерами пересечений подобластей, типами итерируемых краевых условий на смежных границах, количество задействованных подобластей и соответствующих MPI-процессов, число вычислительных потоков на CPU, а также различные счетные параметры для возможного управления скоростью сходимости итераций. Предусмотрены также режимы выборы параметров по умолчанию без вмешательства пользователя.
Заключение
В работе рассмотрены актуальные математические и технологические проблемы решения сверхбольших плохо обусловленных СЛАУ, принципиально влияющие на эффективность параллельных программных реализаций алгоритмов для МВС сложной современной архитектуры. Основное положение заключается в комплексности тесно взаимосвязанных вопросов: обоснование, сходимость и оптимизация крыловских процессов для широкого многообразия типов и свойств матриц; алгебро-геометрические аспекты предобуславлива-ния СЛАУ; алгоритмические и технологические особенности различных видов декомпозиции областей; основные подходы к ускорению крыловских итерационных методов на базе грубосеточной коррекции, дефляции и других приемов улучшения или расширения базисных векторов для подавления ошибки; программные реализации параллельных алгебраических решателей на вычислительных системах с многоуровневой распреденной и общей памятью.
Работа поддержана грантом РФФИ №11-01-00205, а также грантами Президиума РАН №15.9-4 и ОМН РАН №1.3.3-4.
Литература
1. Ильин, В.П. Методы и технологии конечных элементов / В.П. Ильин. — Новосибирск: Изд-во ИВМиМГ СО РАН, 2007.
2. Лебедев, В.И. Вариационные алгоритмы метода разделения области / В.И. Лебедев, В.И. Агошков. — М., Препр. ОВМ РАН; № 54. — 1983.
3. Bramble, J.H. Convergence estimates for product iterative methods with applications to domain decomposition / J. Bramble, J. Pasciak, J. Wang, J. Xu // Math. Comp. — 1991. — Vol. 57, № 195. — P. 1-21.
4. Ильин, В.П. Параллельные методы и технологии декомпозиции областей. / В.П. Ильин // Вестник ЮУрГУ. Серия «Вычислительная математика и информатика». — 2012. — Вып. 1. — №. 46(305). — С. 31-44.
5. Saad, Y. Iterative Methods for Sparse Linear Systems, Second Edition / Y. Saad. — SIAM, 2003.
6. Krylov: библиотека алгоритмов и программ для решения СЛАУ / Д.С. Бутюгин,
В.П. Ильин, Е.А. Ицкович и др. // Современные проблемы математического модели-
рования. Математическое моделирование, численные методы и комплексы программ. Сборник трудов Всероссийских научных молодёжных школ. Ростов-на-Дону: Изд-во Южного федерального университета, 2009. — С. 110-128.
Т. Ильин, В.П. Проблемы высокопроизводительных технологий решения больших разреженных СЛАУ j В.П. Ильин jj Вычислительные методы и программирование. — М., МГУ. — 2009. — Т. 10, № 1. — C. 141-14Т.
В. Бутюгин, Д.С. Методы параллельного решения СЛАУ на системах с распределенной памятью в библиотеке Krylov j Д.С. Бутюгин, В.П. Ильин, Д.В. Перевозкин jj Вестник ЮУрГУ. Серия «Вычислительная математика и информатика». — 2012. — Т. 4Т, № 30б.
— С. 5-19.
9. Intel Math Kernel Library from Intel. URL: http://software.intel.com/sites/ products/documentation/doclib/mkl_sa/11/mklman/index.htm (дата обращения:
15.02.2013).
10. Ильин, В.П. Параллельные методы декомпозиции в пространствах следов / В.П. Ильин, Д.В. Кныш jj Вычислительные методы и программирование. — 2011. — Т. 12, № 1. —
С. 100-109.
11. Brezina, M. An improved convergence analysis of smoothed aggregation algebraic multigrid j M. Brezina, P. Vanek, P.S. Vassilevsky jj Numer. Linear Algebra Appl. — 2012. — Vol. 19.
— P. 441-4б9.
12. Farhat, C. FETI-DP: A dual-primal unified FETI method. Part I: A faster alternative to the two-level FETI method j C. Farhat, M. Lesoinne, P. LeTollei, K. Pierson, D. Rixen jj Int. J. Numer. Math. Engrg. — 2001. — Vol. 50. — P. 1523-1544.
13. Кластер НКС-З0Т: URL: http://www2.sscc.ru/HKC-30T/HKC-30T.htm (дата обращения:
15.02.2013).
14. Message Passing Interface at Open Directory Project: URL: http://www.dmoz.
org/Computers/Parallel_Computing/Programming/Libraries/MPI/ (дата обращения:
15.02.2013).
15. Малышкин, В.Э. Параллельное программирование мультикомпьютеров / В.Э. Малыш-кин, В.Д. Корнеев jj Новосибирск: Изд. НГТУ, 200б.
16. CUDA Tools & Ecosystem: URL: http://developer.nvidia.com/cuda-tools-ecosystem
(дата обращения: 15.02.2013).
1Т. Bell, N. CUSP: Generic parallel algorithms for sparse matrix and graph computations j N. Bell, M. Garland jj URL: http://cusp-library.googlecode.com (дата обращения:
15.02.2013).
1В. Karypis, G. A fast and high quality multilevel scheme for partitioning irregular graphs / G. Karypis, V. Kumar jj SIAM J. Sci. Comp. — 1999. — Vol. 20, № 1. — P. 359-392.
19. Hypre: URL: http://acts.nersc.gov/hypre/ (дата обращения: 15.02.2013).
20. PETSc: Home Page: URL: http://www.mcs.anl.gov/petsc/ (дата обращения:
15.02.2013).
21. Yousef Saad — SOFTWARE: URL: http://www-users.cs.umn.edu/~saad/software/ (дата обращения: 15.02.2013).
Ильин Валерий Павлович, г.н.с. ИВМиМГ СО РАН, д.ф.-м.н., профессор НГУ, e-mail: [email protected]
ON THE QUESTIONS OF PARALLELIZED KRYLOV’S ITERATIVE METHODS
V.P. Il’in, Institute of Computational Mathematics and Mathematical Geophysics of
Siberian Branch of the Russian Academy of Sciences (Novosibirsk, Russian Federation)
Mathematical questions of various computational technologies of parallelelized iterative processes of Krylov’s type for solving large sparse symmetric and non-symmetric SLAEs, obtained in grid approximations of multi-dimensional boundary value problems for PDEs, are considered. Example are presented by finite approximations in gas-hydrodynamical applications, where five unknowns in each node are defined and corresponding SLAEs have small-block structure. The base of used algorithms is flexible generalized minimal residual, FGMRES, method with dynamical preconditioners of additive type, which presents an upper level of two-step iterarive Swartz algorithm.
High performance of algebraic solvers is provided by using different approaches: domain decompositions of various topologies, boundary conditions and sizes of subdomain overlapping, coarse grid correction, deflation and aggregation, and incomplete factorizations of matrices. The unified formulations of using algorithms as well as the questions of computational efficiency and scalable parallelization at the geterogenous supercomputers are described. The examples of technical requirements for peculiarities of program implementation of the libraries of parallel algorithms for solving systems of linear algebraic equation, are presented.
Keywords: iterative methods, Krylov subspaces, preconditioned matrices, domain
decomposition, parallel algorithms, program and computational technologies.
References
1. Il’in V.P. Metody i tekhnologii konechnykh elementov [Finite element methods and technologies]. Novosibirsk, NSC Publ., 2007.
2. Lebedev V.I., Agoshkov V.I. Variazionnyi method dekompozitsii oblastei [Variational domain decomposition method]. Moscow, DCMRAS, preprint No 54, 1984.
3. Bramble J.H., Pasciak J.E., Wang J., Xu J. Convergence estimates for product iterative methods with applications to domain decomposition. Math. Comp. 1991. Vol. 57, No 195. P. 1-21.
4. Il’in V.P. Parallelnye methody i tehnologii dekompozitsii oblastei [Parallel domain decomposition methods and technologies]. Vestnik YUURGU. Seriya “Vychislitelnaya matematika i informatika”[Bulletin of South Ural State University. Seriers: Computational Mathematics and Software Engineering]. 2012. Vol. 46, No 305. P. 31-44.
5. Saad Y. Iterative Methods for Sparse Linear Systems, Second Edition. SIAM, 2003.
6. Butyugin D.S., Il’in V.P., Itskovich E.A. et all. Krylov: biblioteka algoritmov i program dlya teshenia SLAU.[Krylov: the library of algorithms and programs for solving SLAEs]. Modern problems of math. modeling. -Rostov-Don, YUFU Publ., 2009. P. 110-128.
7. Il’in V.P. Problemy vysokoproizvoditelnyh tehnologiy reshenia bolshih redkih SLAU [Problems of high performance technologies of solving large sparse SLAEs] Vychislitelnye
metody i programmirovanie [Computational methods and programming]. 2009. Vol. 10, No 1. P. 141-14Т.
8. Butyugin D.S., Il’in V.P., Perevozkin D.V. Metodu parallelnogo resheniya SLAU na systemah s raspredelennoi pamyatyu [Methods of parallel solving SLAEs on the systems with distributed memory]. Vestnik YUURGU. Seriya “Vychislitelnaya matematika i informatika”[Bulletin of South Ural State University. Seriers: Computational Mathematics and Software Engineering]. 2012. Vol. 4Т, No 30б. P. 5-19.
9. Intel (R) Math Kernel Library from Intel. URL: http://software.intel.com/en-us/
intel-mkl (accessed 12 February 2013).
10. Il’in V.P., Knysh L.V. Parallelnye metody dekompozitsii v prostranstve sledov [Parallel domain decomposition methods in trance subspaces]. Computational methods and programming. 2011. Vol. 12, No 1. P. 100-109.
11. Brezina M., Vanek P., Vassilevsky P.S. An improved convergence analysis of smoothed aggregation algebraic multigrid. Numer. Linear Algebra Appl. 2012. Vol. 19. P. 441-4б9.
12. Farhat C., Lesoinne M., LeTollei P., Pierson K., Rixen D. FETI-DP: A dual-primal unified FETI method. Part I: A faster alternative to the two-level FETI method. Int. J. Numer. Math. Engrg. 2001. Vol. 50. P. 1523-1544.
13. Cluster НКС-З0Т: URL: http://www2.sscc.ru/HKC-30T/HKC-30T.htm (accessed 12
February 2013).
14. Message Passing Interface at Open Directory Project: URL: http://www.dmoz.org/
Computers/Parallel_Computing/Programming/Libraries/MPI/ (accessed 12 February 2013).
15. Malyshkin V.E., Korneev V.D. Parallelnoe programmirovanie multikompyuterov [Parallel programming the multi-computers]. Novosibirsk, NSTU Publ., 200б, 310 p.
16. CUDA Tools & Ecosystem: URL: https://developer.nvidia.com/
cuda-tools-ecosystem (accessed 12 February 2013).
1Т. Bell N., Garland M. CUSP: Generic parallel algorithms for sparse matrix and graph computations: URL: http://cusp-library.googlecode.com (accessed 12 February 2013).
1В. Karypis G., Kumar V. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM J. Sci. Comp. 1999. Vol. 20, No 1. P. 359-392.
19. Hypre: URL: http://acts.nersc.gov/hypre/ (accessed 12 February 2013).
20. PETSc: Home Page: URL: http://www.mcs.anl.gov/petsc/ (accessed 12 February 2013).
21. Yousef Saad - SOFTWARE: URL: http://www-users.cs.umn.edu/~saad/software/
(accessed 12 February 2013).
Поступила в редакцию 13 марта 2013 г.
