Параллельные алгоритмы формирования и решения системы дополнения Шура на графических ускорителях Текст научной статьи по специальности «Математика»

_____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА

Том 154, кн. 3 Физико-математические пауки

2012

УДК 004.272.2^519.612:519.63

ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ ФОРМИРОВАНИЯ И РЕШЕНИЯ СИСТЕМЫ ДОПОЛНЕНИЯ ШУРА НА ГРАФИЧЕСКИХ УСКОРИТЕЛЯХ

С.П. Копы,сов, И.М. Кузьмин,, Н.С. Недож.огии, А.К. Новиков

В работе рассмотрен параллельный алгоритм вычисления дополнения Шура. Эффективное применение нескольких графических ускорителей для метода дополнения Шура связано с разделением матриц и определением алгоритмов, которые более эффективно выполняются па центральном процессоре (СРХГ) или графических ускорителях (СРХГ). Представлен алгоритм обращения матрицы через решение матричной системы множеством параллельных потоков. Показано, что формирование матриц дополнения Шура для нескольких подобластей эффективно выполнять па ЄРи, а с ростом числа подобластей па СРи. Для решения интерфейсной системы предложен параллельный алгоритм метода сопряженных градиентов с явным предобуславливателем. позволяющий достигать существенного ускорения вычислений (в 251 раз) на восьми ЄРІІ при разделении исходной системы уравнений па 64 подобласти.

Ключевые слова: дополнение Шура, параллельные вычисления, предобусловлеп-пый метод сопряженных градиентов, графические ускорители.

В работе рассматриваются параллельные алгоритмы, основанные на блочной декомпозиции матриц и так называемом дополнении Шура. Впервые, метод дополнения Шура был математически сформулирован и описан в работах Исаи Шура (I. БсЬиг), а название своё получил после выхода работы [1]. хотя до этого метод и был известен в вычислительной механике как метод подструктур и впервые встречается в работе [2]. В нашей стране он известен как метод суперэлементов [3]. Метод дополнения Шура используется также для обращения больших матриц при помощи разбиения на клетки [4]. Первоначально [2] он рассматривался как метод декомпозиции области для вычислений на системах с ограниченными вычислительными мощностями. В дальнейшем этот метод использовался для построения методов декомпозиции области без перекрытий при разработке параллельных алгоритмов решения задач вычислительной механики. В настоящее время его всё чаще рассматривают как гибридный метод решения систем линейных алгебраических уравнений (СЛАУ) [5. 6]. сочетающий преимущества как прямых, так и итерационных методов и учитывающий различные архитектуры параллельных вычислительных систем. Кроме того, на основе дополнения Шура строятся эффективные предобуславливатели для решения СЛАУ различных видов (см., например. [7]).

Во многих из указанных задач требуется решить блочную (2 х 2) систему уравнений вида

Аннотация

Введение

Дополнение Шура матрицы А11 определяется как

& — А22 — А21Ац1 А12?

если матрица Лц € М”х” - квадратная и невырожденная, А12 € М”хт, А21 € € Мтх”, А22 € Ктхт. Тогда редуцированная система уравнений для нахождения Х2 записывается в виде

5X2 = /2 - А21А/

Неизвестные Х1 находятся из соотношений

Х1 = А11/1 - А12А111Х2.

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

Современные графические ускорители (СР11) позволяют для ряда задач, в том числе при работе с векторами и матрицами, получать ускорение вычислений в десятки и сотни раз по сравнению с центральным процессором (СР11). Тысячи потоков (нитей) СРи могут эффективно выполнять одновременно большое число простых арифметических операций, что характерно для мультипликативных и аддитивных операций с векторами и матрицами. Вместе с тем последовательные операции и ветвления, характерные, например, при разложении матриц на треугольные множители, выполняются на СР11 медленнее, чем ядрами центрального процессора.

Таким образом, эффективное применение СР11, тем более нескольких, для метода дополнения Шура связано с декомпозицией матриц и определением алгоритмов. которые более эффективно выполняются на СР11 или на графическом ускорителе.

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

1. Ресурсная эффективность метода дополнения Шура

Рассмотрим процесс построения дополнения Шура как один из вариантов методов декомпозиции области в виде алгоритма метода подструктур [2]. Для обеспечения независимости вычислений в отдельных подобластях и последующего их взаимодействия все узлы области делятся на два множества: внешние и внутренние. Неизвестные перемещения области рассматриваются в виде суперпозиции двух составляющих. Первая составляющая перемещения, вызванные внешними силами при закреплении границ в подобластях. Перемещения каждой подобласти определяются из уравнений, включающих неизвестные, связанные только с данной подобластью. Вторая составляющая перемещения, вызванные смещениями границ подобласти с исключенными внутренними узлами.

Пусть область П разбита на пп непересекающихся подобластей:

пп

П = П^ П2и Ппп, оде П<р| Пй = 0, Гв = и дП \ дП. (2)

£=1

Разделение на подобласти наследуется от процесса разделения дуального графа

пп

расчетной сетки С(У, Е) = |^| С£(Т£, Е£), здесь множество V вершин графа - это

£=1

Е

жество смежных конечных элементов. Множество конечных элементов, образующих подобласть, - это V С V. Далее полагается, что все подграфы С£(Т£,Е£) связные, в противном случае система уравнений для П£ подобласти распадается на несвязанные между собой системы уравнений.

Узлы расчетной сетки образуют множество V и условно разделяются на внешние Ув* (принадлежат границе области) и внутренние Т7^ (связанные с узлами подобласти) узлы сетки, соответствующей подграфу С£(Т£,Е£). Из множества внешних узлов выделяются интерфейсные С ув*, связанные с узлами из других подобластей. В соответствии с этими обозначениями и введённым упорядочиванием система уравнений примет вид

Л *-1 Л >-ч СО >-ч Л|в\ Л2в ии *“ч СО>~ч 1—*■ 1-ЧкЦкч //

^Лв/ Л2в1 ■ Л Л (л 31 г Лпп Л1в Лвв/ ипп ив /пп /в

где верхний индекс г определяет принадлежность подобластн П £, а нижние индексы I, В относятся к внутренним и граничным степеням свободы. Здесь невы-

пп

рожденная матрица Лц = ЛгП соответствует А11 из системы уравнений (1),

£=1

пп

матрица Авв - блок у Л22, блок матри цы Л1в = Л}в - блок у Л12. Отме-

£=1

тим, что блок матрицы Лц положительно определен и симметричен, так же как и полная матрица в (3) для рассматриваемого класса задач.

Система для интерфейсных узлов есть

Явв ив = /в. (4)

Пп / -1 \

Здесь Бвв =^Г^ [Лгвв — Лгв1 Л\1 Л}в ) - матрица дополнения Шура для подобла-

г

- Пп / _1 \

сти г, вектор /в = (/в — Л1в1 Л111 /А - вектор правых частей.

г

Для решения задач используется усовершенствованный алгоритм метода подструктур. В его основе лежит использование свойств невырожденности и положительной определенности матриц подобластей. Для таких матриц существует разложение Холесского Лц = Ьц ЬЦ, где Ь - нижняя треугольная матрица с положительными диагональными элементами. Использование разложения Холесского значительно сокращает вычислительные затраты и объем используемой оперативной памяти вычислительных машин.

Ниже рассмотрена реализация последовательного алгоритма вычисления дополнения Шура (пп > пр и число процессоров пр = 1) и приведены вычислительные затраты, связанные с каждым шагом алгоритма (для каждого шага

алгоритма показано число необходимых операций с учетом симметрии матриц, причем операции сложения и умножения принимаются за одну).

Алгоритм 1. Последовательный вариант дополнения Шура

1: Выполним разложение Холесского матрицы Ajj для соответствующей подобласти

(индекс i опущен)

Лц = Ьц Lfj . (n 3/6)

'2: Вычисляем вспомогательные переменные

AjB = Л//Л1В ■ (пЕ • п|)

3: Сформируем матрицы граничной жесткости подобластей

Sbb = Лвв — AbjAjb . ((nj • пЕ)/2)

4: Формируем вектор правой части

/в = /в — AbjЛ-/ fj . (nj • ПЕ )

о: Собираем и решаем систему уравнений

Sbbub = /в ■ (k(M-1Sb Е) < C(1+log(H/h)))

6: Определяем неизвестные для внутренних узлов

uj = Л-/ fj — Ajb ub . (»?)■

Приняты следующие обозначения: n/ = m • |Vffc I и пв = m • |Увк | — число внутренних и граничных степеней свободы соответственно, m - число степеней свободы в узле сетки, h - шаг сетки, H - размер подобласти, M - предобуславливатель для дополнения Шура. На пятом шаге Алгоритма 1 приведена оценка обусловленности матрицы SBB, а не вычислительная сложность, которая зависит от выбора метода решения системы уравнений. В настоящей работе в качестве метода решения СЛАУ использовался метод сопряженных градиентов с различными предобу-славливателями. Матрицы Sbb j A// положительно определены и симметричны. Отметим, что порядок матрицы Sbb значительно меньше, чем исходной матрицы, но достаточно большой. Поэтому для решения системы SBBuB = fB применяются итерационные методы решения и компактные схемы хранения матриц.

Для обеспечения эффективной работы с матрицами используются различные схемы хранения. Матрицы Sbb > A// являются симметричными, поэтому возможно хранение только её части (верхней или нижней треугольной и главной диагоналей), оценка ресурсоёмкости схем хранения производится в зависимости от способа хранения симметричных матриц. Простейший вариант хранить матрицу, не исключая нули. В этом случае для хранения используется массив по числу элементов матрицы N2. Этот способ является наиболее затратным с точки зрения требуемой памяти, что особенно существенно для больших задач.

Формат хранения DCSR, используемый при вычислениях в данной работе, представляет собой массив, состоящий из списка упакованных строк, реализован в системе конечноэлементного анализа FEStudio [8, 9]. Каждая строка матрицы представляет из себя структуру, состоящую из двух массивов. В первом массиве хранятся значения ненулевых элементов, во втором столбцовые индексы этих элементов. Размер каждого из массивов для строки i равен числу ненулевых элементов Nnz и меняется динамически по мере необходимости.

Предложенный формат DCSR является разновидностью распространенного формата хранения разреженных матриц формат CSR (Compressed Sparse Row Storage Format). Для матрицы A при использовании формата CSR выделяются три одномерных массива, в которых хранятся ненулевые значения {aj | aj = 0}, 1 ^ ij ^ Nnz, их столбцовые индексы {j | aj =0}, 1 ^ i,j ^ Nnz, и номера элементов aji, 1 ^ i ^ N + 1,в двух первых массивах.

Оценим рссурсоемкость предложенной схемы хранения матриц, для чего сравним три формата хранения матриц (ленточный. БСБК и СБИ) по следующим параметрам: требуемая память, алгоритмическая сложность доступа к элементу и его внесения. Сравнение показывает, что алгоритмические сложности доступа к элементу 0(2Жд + 1) и его внесения для ленточного формата — 0((2Жд + 1)Ж). для БСБК - ©(N«2*) и ©(N«2*) соответственно, а для формата СБИ - ©(N«2*)

и ©(N«2), здесь N«.2* = тах \Nnzi}. Необходимый объем памяти для ленточ-

i=1,...,N

N

ного формата - 8(2Жд + 1^, для формата БСБК - + 4), для формата

І= 1

СБИ - 12Nnz + 4^ + 1) байт. Для для формата БСБК полагается, что при храпении вещественной величины используется 8 байт памяти и 4 байта для целого числа. В симметричном случае величины 2Жд + 1 в оценках сложности алгоритма (и приведенной далее частоты доступа) заменяются на N0 + 1, а N«2 и N«2.* относятся к элементам треугольной матрицы.

Отметим, что формат БСБК более удобен, чем СБИ, при выполнении треугольного разложения матриц Лц, когда в строках появляются новые ненулевые элементы. В этом случае изменение в одной нз строк не требует смещения всех последующих элементов в массивах, как в формате СБИ. Поэтому процедура добавления нового элемента имеет меньшую алгоритмическую сложность ©(N«2*), чем в формате СБИ - ©(N«2). Кроме того, среди представленных форматов объем памяти,

N

необходимый в случае формата БСБК, равный У^(12Nnzi + 4) байт, является

І=1

минимальным.

Преимущества ленточного формата по сложности доступа и добавления элемента нивелируются увеличением необходимого объема памяти 8(2Жд + 1^ и частоты доступа к элементам матрицы ©((2Жд + 1^) вместо ©(N«2) для форматов БСБК и СБИ. Здесь Np = тах^=1 тах^О' — * | а^ = 0} - так называемая полуширина ленты, зависящая от нумерации неизвестных и уравнений, N - размерность матрицы, N«2 - число ненулевых элементов в матрице, Nnzi - число ненулевых элементов в >-й строке матрицы.

Как показали эксперименты [10], схема хранения оказывает существенное влияние на время вычислений. Отметим, что для ленточных матриц время формирования дополнения Шура составляет порядка 80% от общего времени, затраченного на решение. Переход на формат БСБК для матриц Лвв, Л1в, Лц для одной и той же задачи дал сокращение затрат в четыре раза. Таким образом, в последовательном варианте наиболее эффективным, с точки зрения ресурсоемкости и алгоритмической сложности, представляется использование метода дополнения Шура с разложением Холесского матрицы Лц и храпением матриц в сжатом формате.

Обработка строк матрицы, хранящейся в формате БСБЕ на графическом ускорителе (СР11), потребует для каждой строки: выделения памяти на СР11, копирования строки, далее в ядре выполнения вычислений над строкой и возвращения результата в оперативную память или кэш СР11. Таким образом, вместо выделения памяти и копирования двух массивов размера N«2 потребуется 2N выделений памяти и копирований массивов размера Nnzi, поэтому для вычислений на СРи матрица из формата Б СБИ переводится в формат СБИ.

2. Параллельная эффективность метода дополнения Шура

В методе декомпозиции области можно выделить два направления распараллеливания процесса решения задачи. Первое связано с методами явного деления на подобласти, когда параллельные вычисления осуществляются на высоком уровне

(уровне области), то есть число ветвей параллельного алгоритма равно числу подобластей. Как показывает практика, такой подход дает значительное ускорение, если каждый процессор решает в своей подобласти свою подзадачу. Второе направление (уровень подобласти) неявные методы подобластей. Внутри каждой подобласти различные задачи могут быть решены наиболее эффективным спосо-

Рассмотрим. исходя из обозначенных вариантов распараллеливания. Алгоритм 1. Как видно. Шаги 1 4, 6 выполняются для каждой подобласти независимо, что говорит о естественном параллелизме, который обеспечивает первый уровень распараллеливания на уровне области. Распараллеливание Шага 5 зависит от выбранного метода решения системы (4).

Наиболее трудоёмкими являются шаги формирования (Шаги 1 4) н решения системы дополнения Шура (Шаг 5). Поэтому распределение вычислительных затрат на этих шагах по нескольким процессорам (центральный процессор, графический ускоритель) является решающим фактором для параллельного ускорения вычислений и определяет эффективность распараллеливания Алгоритма 1.

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

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

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

Так. сбалансированное распределение вычислительной нагрузки при выполнении Шагов 1~4 зависит не от общего числа узлов |у | в подобластях Qi, а от числа внутренних |У1 | и внешних |Ув* | узлов в подобластях.

Вместе с тем сетки для трёхмерных областей с развитой поверхностью (пружины. тонкие пластины и оболочки) содержат большое число конечных элементов. у которых все узлы принадлежат границе расчётной области. Разделение дуальных графов таких сеток и выполнение условия |у| « |у |5 * = 3, приводит к подобластям ^, та которых |у4 | = 0. Наглядным примером является задача моделирования напряженно-деформированного состояния пружины, рассматриваемая в настоящей работе. Геометрия пружины аппроксимируется неструктурированной расчетной сеткой (см. рис. 1, о) с ячейками в виде тетраэдров, число ячеек |У| = = 174264, число узлов |у| = 40743. Для получения подобластей па которых |У1 | > 0, дуальный граф С(У, Е, Ш) полагался взвешенным с множеством весов Ш = {и^}. Если хотя бы один из узлов конечного элемента не лежит на дП, то вес соответствующей вершины графа ик = 3, в противном случае ик = =1

(«п = 16, 32,..., 1024). Отметим, что число подобластей увеличилось в 64 раза,

а) б) о)

Рис. 1. Сетка: а) исходная: 6) разделенная па 16 подобластей: а) разделенная па 1024 подобласти с помощью двухуровневого разделения

а размер матрицы £вв - только в 1.44 раза (Ж = 66030 в случае «п = 16 и N = 95523 для «п = 1024). При этом число ненулевых элементов матрицы дополнения Шура уменьшилось в 18 раз (с N«2 « 2.7 • 108 в случае 16 подобластей, до N«2 « 1.5 • 107 - для 1024 подобластей), как следствие, уменьшилась её заполненность с 6.11% до 0.16%. В среднем примерно каждый шестнадцатый элемент в строке ЙВд является ненулевым (4041 из 66030) при «п = 16 (см. рис. 1, б) и примерно каждый шестисотый (154 из 95523) при разделении на 1024 подобласти (см. рис. 1, е).

Важно отметить, что размер системы для дополнения Шура меньше размера исходной конечноэлементной системы (в рассмотренных случаях в 1.4 2.0 раза). При этом заполненность матрицы на один-два порядка больше, чем запол-

ненность глобальной матрицы жесткости (0.03%).

3. Формирование матрицы дополнения Шура

Одной из самых затратных операций формирования дополнения Шура является обращение матрицы Ац. Обычно на практике для нахождения обратной матрицы используются плохо распараллеливаемые прямые методы, например метод обращения на основе -разложения.

Рассматриваемый в работе алгоритм вычисления обратной матрицы состоит из решений матричной системы вида АцX = Е, где Е - единичная матрица «/ х «/. Система эффективно решается на СР11 иредобусловленным алгоритмом сопряженных градиентов (см. следующий параграф). Если в правой части системы вместо матрицы Е брать А1в, то её решением будет матрица А1в = А-/А/В. Такое представление позволяет заменить операции обращения матрицы и матричного произведения па решение «в систем, каждая из которых решается независимо, что позволяет использовать одновременно несколько СРи. Дополнение Шура, или матрица граничных жесткостей, вычисляется по соотношениям (4), где Авв € € М”в х”в, Ап € М”1 хиг, Ав/ € М”в хиг, А/в € М”1 х”в .

Пусть пВ ^ количество столбцов матрицы Авв пересылаемых па *-й СР11,

АВВ € М”в х”в - матрица, состоящая из столбцов матрицы Авв с номерами от І І+1

13 до И пВ , где * - номер СРи. Каждый графический ускоритель решает пгв

5=0 5=0

Табл. 1

Время формирования дополнения Шура, мин

пп сри і єрі: 2 єрі: 4 єрі: 6 єрі: 8 єрі:

16 26:23.6 31:52.6 19:25.5 13:09.9 10:57.6 09:49.4

32 08:00.6 19:33.9 11:01.8 06:41.7 05:14.0 04:26.5

64 02:25.1 11:18.8 06:18.2 03:44.8 02:54.7 02:29.0

128 03:57.2 02:25.8 01:58.5 01:42.5

256 03:14.0 02:01.0 01:37.7 01:26.0

512 02:48.3 01:45.7 01:26.0 01:15.4

1024 00:18.7 03:53.4 02:24.0 01:32.2 01:17.5 01:08.0

систем Ацай = а^в , здесь а^в - к-й столбец матрицы А/в, к €

А/в = {а1, а2 ... а”в }. В реализации подхода независимого решения систем уравнений на нескольких графических ускорителях используется технология ОрепМР. Для этого создаются несколько нитей (число которых равно количеству доступных устройств СР11). определяется номер нити и каждой назначается графический ускоритель с тем же номером.

Ниже представлен Алгоритм 2, реализующий этот подход.

Алгоритм 2. Параллельный алгоритм формирования дополнения Шура на каждой подобласти П (индекс і у матриц опущен)

1: Решаем систему ЛііЛ'ІВ = Лів; {матрицы хранятся на СРІІ в СвІІ-формате, решение находится по столбцам}.

2: Явв = Лвв — ЛвіЛів ; { Явв и результат произведения матриц - построчно, Лвв -СвІІ-формате}

3: fв = /в — Лвіх; {х _ решение системы Лііх = /і}.

4: Формируем и решаем: Явв ив = /в; {Явв формируется в БСБЕ на СРи, а затем копируется в СвІІ па СРП }.

5: Определяем иі = х — Лів ив ; {ха Лів вычислены ранее и хранятся на СРІІ}.

На каждом ЄРІІ после решения систем остаётся матрица (А1В) Є М” х”в , со-

І і+1

стоящая из столбцов матрицы А1В с номерами от Е П-В до ЕПВ . Это позволяет

5=0 5=0

без дополнительных коммуникаций выполнить оставшиеся операции (произведение и разность матриц, см. соотношения (4)) для каждой матрицы (А1В)і независимо на нескольких ЄРІІ, которые используются при вычислении матриц £Вв • Далее формируется глобальная матрица дополнения Шура £ВВ (Шаг 4 в Алгоритме 2), состоящая из локальных £Вв , принадлежащих і-й подобласти.

Локальные матрицы дополнения Шура ^ВВ хранятся в несжатом формате. Матрица £ВВ после формирования из локальных £ВВ преобразуется из несжатого к тому формату, в котором будет наиболее удобно с ней работать на СРІІ, например к СБТІ-формату.

В табл. 1 приведены временные затраты по параллельному формированию матриц дополнения Шура в зависимости от числа подобластей и СРІІ. Во второй колонке приведены данные для последовательного алгоритма, выполняемого на центральном процессоре (Алгоритм 1). Ускорение параллельного алгоритма (Алгоритм 2), реализованного на СРІІ, по отношению к СРІІ, определим как в(пр)сри = = £сриД(пР)ори > гДе ^СРи _ время выполнения последовательного алгоритма, реа-

ЕПВ

І+1

Е5

лизованного на СРІІ, а і(пр)сри - время выполнения параллельной реализации на пр ЄРи. Аналогичным образом введём ускорение в(пр)сри = і(1)вриД(пр)сри-В рассмотренных вариантах максимальное ускорение составляет в(8)ори = 2.7 и достигается при числе подобластей п^ = 16, при этом в каждой подобласти находится в среднем примерно по 11000 ячеек сетки. Формирование матрицы дополнения Шура на двух СРІІ приводит к ускорению в(2)ори > 1.5 в зависимости от числа подобластей. Использование восьми графических ускорителей даёт наименьшее время выполнения для реализации на СРІІ, но для задач с небольшим числом ячеек в каждой подобласти (< 2500), СРІІ-реализация оказывается эффективнее. В этом случае при решении большого числа систем уравнений малой размерности для каждой подобласти требуется время на копирование данных между СРІІ и СРІІ, которое становится большим, чем время решения систем. Для одного СРІІ (при использовании нескольких) затраты на решение сокращаются, но не покрывают затрат на инициализацию и копирование.

Полученные данные позволяют применять рассматриваемый алгоритм формирования для обеспечения более сбалансированного использования СРІІ и СРІІ в гибридных вычислительных системах.

Матрица дополнения Шура € М№ х№ (далее S = [в^ ]) имеет порядок

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

В большинстве работ, посвящённых реализации итерационных методов на СР11, рассматривается предобуславливатель Якоби или его блочный аналог. Оптимальным выбором для вычислений на СР11 нам представляются предобуславлнвателн, для которых предполагается, что известна аппроксимация обратной матрицы системы [11]. В этом случае дополнительные операции, вызванные переходом к пре-добусловлеиной системе, сводятся к матрично-векторному произведению =

В настоящей работе рассматриваются методы сопряжённых градиентов с диагональным предобуславливателем (В1АС), предобу сдавливанием по методу симметричной последовательной верхней релаксации (ББОК) и предобуславливателем, построенным масштабированием системы (БЕР). К сожалению, в наших экспериментах предобуславлнвателн, основанные на аппроксимации обратной матрицы с использованием дополнения Шура, не показали сравнимую эффективность вычислений на рассматриваемых матрицах.

Диагональный предобуславливатель имеет вид:

Определим предобуславливатель ББОК следующим образом. Пусть матрица системы представима в виде Б = Ь + О + Ьт, тогда

4. Решение интерфейсной системы уравнений на СРП

Мгй+і.

если і = л;

.

или

Алгоритм 3. Алгоритм метода сопряжённых градиентов с прсдобуславливатслсм

1: S, М Є lJVxJV {М формируется па GPU, матрицы хранятся в CSR-формате}

2: u, т, p, q, z Є RN {Вектора хранятся в памяти GPU, копии на CPU нет}

3: r0 ^ f {копирование векторов осуществляется с помощью cublasDcopy}

4: u0 ^ 0 {инициализация выполняется па GPU} о: Zo <— Мго {выполняется па GPU}

fi: p0 ^ z0 {копирование векторов осуществляется с помощью cublasDcopy}

7: p0 ^ (т0, z0) {здесь и далее (•, •) = 0(P^}

P

8: while 11Ті ||2/1|bII2 > є do

9: qt ^ Spt {выполняется на GPU}

10: at ^ (Ti,zi)/(qi,pt) {вычисляется с помощью функции cublasDdot}

11: ut+i ^ ut + atpt {операция выполняется с помощью функции cublasDasxpy}

12: Гі+і <— /’г — {операция выполняется с помощью функции cublasDasxpy}

13: Zi+1 <— Мгі+і {выполняется па GPU}

14: pt+i ^ (rt+i, zt+i) {вычисляется с помощью функции cublasDdot}

15: et+i ^ pi+i/pi

16: pt+i ^ zt+i + ві+ipt {последовательное использование функций cublasDscal и

cublasDasxpy}

17: end while

где 0 < w < 2, D = (1/w)D. Приближенно обратную матрицу вычислим, ограничиваясь первым членом в разложении

К = К-1 ^ V2^D-1/2{I - LD-1). (5)

Тогда предобуславливатель SSOR примет вид

м = ктк.

Если положить w = 1 в (5), получим предобуславливатель вида

K = D-1/2(I - LD-1)

и, соответственно,

М = D^2KKTD^2.

При построении иредобуславливателя DIP масштабируем исходную систему:

S = D-1/2SD-1/2; f = D-1/2f; и = D1/2u;

тогда предобуславливатель примет вид

M = (I - LT)(I - L).

Решение интерфейсной системы методом сопряжённых градиентов распараллелено с помощью технологии CUDA для вычисления на GPU. Все вспомогательные массивы, в частности г, p g, z, а также матрица системы, предобуславливатель, вектор правых частей и вектор решения хранятся в памяти графического ускорителя (см. Алгоритм 3). После завершения работы метода сопряжённых градиентов, массив и, в котором хранится приближение вектора решения, копируется в память CPU.

Для реализации операций суммы, скалярного произведения, копирования векторов и умножения вектора на скаляр использовались функции, предоставляемые

библиотекой СиВЬАБ. При выполнении матрично-векторного произведения, вектор хранится в текстурной памяти, которая кэшируется, что даёт более быстрый доступ н уменьшает временные затраты. Для вычисления каждой координаты вектора результата используется от 2 до 32 потоков в зависимости от разреженности матрицы. Вычисление предобуславливателя выполняется либо на СР11, либо на СР11 в зависимости от типа предобуславливателя.

Предобуславливатель ББОК вычисляется на центральном процессоре. Разреженность матрицы К равна разреженности матрицы Б, и поэтому выгоднее её также хранить в компактном формате. Теоретически метод сопряжённых градиентов решает систему не более чем за N итераций. Пусть система решается за N

—т — '

итераций, значит, при вычислении К К требуется N матрично-векторных произведений плюс N матрично-векторных произведений = Мг^\ для решения

системы методом сопряжённых градиентов. Если не вычислять М явно, а произведение = Мгк+1 заменить на два следующих друг за другом матричновекторных произведения = Кгк+1 и ^6+1 = К 2^+1, то нам потребуется

всего 2N матрично-векторных произведений, что меньше, чем N + N — при нахождении матрицы М в явном виде. Поэтому при реализации метода сопряжённых градиентов с ЗБОК-иредобуславливателем произведение ^+1 = Мгк+1 заменено иагк+1 = К Кгк+1.

Масштабирование системы и вычисление предобуславливателя Б ЕР выполняется на СРи. Как и в случае с ББОК, разреженности матриц совпадают, и они также хранятся в разреженном формате. Произведение ^+1 = Мгк+1 также заменяем на два последовательных матрично-векторных произведения у = (I — Ь)г^+1 и 2^+1 = (I — Ьт)у, выполняемых на СР!1.

Применение предобуславливателя ББОК сокращает число итераций при решении системы в полтора раза (оптимальный параметр релаксации для данной задачи ш = 0.5), что позволяет покрыть затраты на его создание и выполнение дополнительного матричио-векториого произведения. В результате этого нет выигрыша по времени вычислений в сравнении с диагональным предобуславливателем. Аналогичные результаты получены для систем разрывного метода Галёркина [11].

Однако отличие в свойствах матрицы Б и матриц систем разрывного метода Галёркина позволило значительно уменьшить число итераций при использовании предобуславливателя Б1Р, которое стало примерно равным их числу в случае использования предобуславливателя ББОК.

Для решения интерфейсной системы (4) реализован блочный алгоритм сопряженных градиентов, использующий пр графических ускорителей. При разделении на блоки , к = 1,..., пр, матрице Б системы (4) соответствовал граф, имеющий число вершин, равное размерности Б. Каждой вершине графа матрицы Б, на основе вычисленного разделения многоуровневым алгоритмом [12], ставится в соответствие номер графического ускорителя (далее в терминах раскрашивания графов «цвет»). Вершины раскрашенного графа подразделяются на внутренние и граничные (связанные хотя бы с одной вершиной, имеющей другой «цвет»).

Па основе полученного разделения в каждом блоке выделялось несколько

гч\ЪкАк] «

матриц: Б к — матрица, элементы которой связывают внутренние вершины в блоке; ’Ьк], Б^к,1к] — матрицы, связывающие внутренние вершины с граничными; Б\Ьк,Ьт] - матрица, связывающая граничные вершины к-го блока с граничными вершинами т-го блока. Здесь к = т и к,т € [1,пр], где пр - число блоков.

Табл. 2

Время решения интерфейсной системы. ч:мип:с

Шг сри і єрі: 2 єрі: 4 єрі: 6 єрі: 8 єрі:

16 > 8:00:00 0:15:12 0:07:01 0:03:25 0:04:41 0:04:24

32 > 8:00:00 0:16:23 0:10:30 0:07:10 0:05:17 0:04:34

64 5:01:07 0:04:47 0:02:54 0:01:38 0:01:22 0:01:12

128 0:02:07 0:01:18 0:01:06 0:01:00

256 0:01:46 0:01:10 0:01:01 0:00:56

512 0:01:25 0:01:03 0:00:56 0:00:53

1024 1:48:22 0:01:09 0:01:16 0:00:59 0:00:54 0:00:52

Запишем матрицу S в виде

0 0

0

V 0

[*і,Ьі.

[Ьі,Ьі]

[Ьь,Ьі]

0

0

[ік ,гк] к

[Ьк,ік_

[Ьі,Ьь]

'[ік ,Ьк_ к

■[Ьь,Ьь]

к

0

[Ьі.Ьпр] 1

0

[Ьк,Ь„

0

[Ьп„ ,Ьі

0

[Ьп„ ,Ьк]

С,[*п

[Ь„

£,[іпр ?Ьпр ] 5Пр

5 [Ьпр ,Ь„р] 5"р /

При выполнении матрично-векторного произведения матрицы 5 на вектор

рт = ^р!, р 1,..., рк, Рік, • • •, РПр, Рпр) на каждом СРи вычисляются два подвек-тора

т^Пр т= 1

V •

т

уг = 5 [^'^]рг і 5 [ік'Ьк]РЬ ук = 5 к Р к + 5 к Р к,

(6)

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

Решение СЛАУ методом сопряженных градиентов требует меньше временных затрат в случае алгоритма, использующего СРІІ (см. табл. 2), при этом ускорение составляет в(1)сри = 72 при разделении та 16 подобластей, и в(1)сри = 94 при разделении на 1024 подобласти. Использование нескольких СРІІ даёт максимальные ускорения в(8)сри = 251 для п^ = 64 и в(8)сри = 3.5 для п^ = 16. С увеличением числа подобластей количество ненулевых элементов в полученной матрице дополнения Шура уменьшается, это приводит к тому, что эффективность использования нескольких СРІІ для решения интерфейсной системы снижается. Так, например, при разделении на 1024 подобласти в(8)ори = 1.3. Шаг 6 в Алгоритме 1, отвечающий за нахождения решений на внутренних узлах, также выполняется на ЄРи. Произведение А—1 / уже вычнсленио на Шаге 4, матрица получена на Шаге 2. Поэтому необходимо сделать лишь две операции: матричновекторное произведение и разность векторов, которые выполняются на СРІІ.

В ходе проведення вычислительных экспериментов минимальные значення суммарного времени формирования и решения системы для дополнения Шура (4) получены при п^ = 1024 (см. табл. 1 и 2). При выполнении этих шагов только

1

центральным процессором потребовался 1 ч 48 мин. В случае использования только одного GPU для формирования и решения СЛАУ затраты сократились в 22 раза. Минимальное время вычислений на графических ускорителях получено на восьми GPU: £(8)gpu = 2 мин. Формирование системы (4) на центральном процессоре и решение на одном графическом ускорителе выполнено за полторы минуты. Наименьшие суммарные затраты потребовались при формировании системы на CPU и её решении на восьми GPU.

Полученные результаты показывают, что при решении задач с помощью метода. использующего дополнения Шура, оптимальный выбор алгоритма зависит от числа и размера подобластей, на которые разбивается расчётная сетка. Если в одной подобласти находится относительно небольшое число ячеек сетки (< 5000), то для формирования матрицы дополнения Шура эффективнее использовать прямые методы нахождения обратных матриц и. как следствие, задействовать только CPU. При больших размерах подобластей наиболее эффективны итерационные алгоритмы. использующие для вычислений несколько GPU. Решение системы уравнений дополнений Шура эффективнее производить на GPU. В этом случае время решения сокращается в десятки и сотни раз.

Метод дополнения Шура следует применять для обеспечения более сбалансированного использования GPU и CPU на гибридных вычислительных системах при решении СЛАУ и в других приложениях.

Отметим некоторые нз возможных обобщений предложенных алгоритмов. Все представленные в работе алгоритмы обобщаются без изменений на случай, когда подматрица А22 в (1) является нулевой. Системы такого типа возникают, например, при использовании смешанных аппроксимаций для численного решения эллиптических уравнений второго порядка. В случае несимметричной матрицы дополнения Шура при решения системы (4) может применяться метод бнеопря-женных градиентов с матрично-векторными произведениями в виде (6). Для обращения несимметричных подматриц Ац в рамках Алгоритмов 1, 2 разложение Холесского необходимо заменить на LU-факторизацию. Для систем уравнений, в которых блок Ац вырожден, предложенный алгоритм обращения может быть заменён на псевдообращение Мура Пенроуза, остальные реализации матричных операций на CPU или GPU остаются без изменений.

Вычислительные эксперименты выполнены на узлах гибридного кластера «Уран» ИММ УрО РАН, каждый из которых содержит два процессора Intel Хеоп Е5675, восемь графических ускорителей NVIDIA Tesla М2090.

Работа выполнена в рамках программы Президиума РАН X- 18 при поддержке УрО РАН (проект 12-П-1-1005) и РФФИ (проект ЛМ1-01-00275-а).

Summary

S.P. Kopysov, I.M. Kuzmin, N.S. Nedozhogin, A.K. Novikov. Parallel Algorithms for Constructing and Solving the Scliur Complement 011 Graphics Accelerators.

The paper deals with a parallel algorithm for computing the Scliur complement 011 multiple GPU. The implementation of a parallel subdomain is shown at the stage of constructing the Scliur complement matrices. A11 algorithm for matrix inversion is presented by the solution of the matrix system for multiple parallel streams. The realization of the mat.rix-vect.or product by means of the matrix decomposition algorithm is described for a parallel conjugate gradient method proposed for the interface system solution.

Key words: Scliur complement., parallel computing, preconditioned conjugate gradient method, graphics accelerators.

Литература

1. Haynsworth E.V. On t.lie Scliur Complement // Basel Math. Notes. 1968. No 20.

17 p.

2. Przcmicniccki J.S. Theory of Matrix Structural Analysis. N. Y.: McGaw-Hill. 1968. 480 p.

3. Постное В,А, Метод суперэлемептов в расчетах инженерных сооружений. Л.: Судостроение. 1979. 288 с.

4. Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. М.: Физ-

матгиз, 1960. 656 с.

5. Girautl L., Haidar A., Saad Y. Sparse approximations of the Scliur complement for parallel algebraic hybrid solvers in 3D // Numer. Math. 2010. V. 3. P. 276 294.

6. Rajamanickam S., Boman E.G., Hcroux M.A. SliyLU: A Hybrid-Hybrid Solver for Mult.icore Platforms // IEEE 26t.li Int. Parallel and Distributed Processing Symposium (IPDPS), 21 25 May 2012. P. 631 643.

7. Kopueeo В.Г., Encen G., Эффективное предобуславливапие методом декомпозиции

области для р-версии с иерархическим базисом // Изв. вузов. Матем. 1999. Л'! 5.

С. 37 56.

8. Копы,сое С.П., Красноперое И.В., Рынков В.Н. Объектпо-ориептироваппый метод декомпозиции области // Вычислительные методы и программирование. 2003.

Т. 4, Л» 1. С. 176 193.

9. Kupysuv S.P., Krasnopyorov I.V., Novikov А.К., Rychkov V.N. Parallel Distributed Object-Oriented Framework for Domain Decomposition // Domain Decomposition Methods in Science and Engineering. Springer. 2005. V. 40. P. 605 614.

10. Копы,сов С.П. Оптимальное разделение области для параллельного метода подструктур // Сеточные методы для решения краевых задач и приложения. Материалы Пятого Всерос. семинара. Казань: Изд-во Казан, ун-та, 2004. С. 121 124.

11. Копы,сов С.П., Новиков А.К., Сагдеева, Ю.А. Решение систем уравнений метода

Галёркипа с разрывными базисными функциями па графическом ускорителе // Вести. Удмурт, уп-та. Математика. Механика. Компьютерные пауки. 2011. Л'! 3.

С. 137 147

12. Karypis G., Kumar V. Parallel multilevel k-way partitioning scheme for irregular graphs // SIAM Rev. 1999. V. 41, No 2. P. 278 300.

Поступила в редакцию 18.06.12

Копысов Сергей Петрович доктор физико-математических паук, профессор, заведующий лабораторией Института механики УрО РАН, г. Екатеринбург.

E-mail: s.kopysovegmail.com

Кузьмин Игорь Михайлович младший научный сотрудник Института механики УрО РАН, г. Екатеринбург.

E-mail: imkuzmin Qgmail. сот

Недожогин Никита Сергеевич аспирант Института механики Уро РАН, г. Екатеринбург.

E-mail: NcdozhoginMinbox.ru

Новиков Александр Константинович кандидат физико-математических паук, старший паучпый сотрудник Института механики УрО РАН, г. Екатеринбург.

E-mail: sc work0mail.ru