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

Уфа : УГАТУ, 2013

УДК 519.63

Параллельные алгоритмы метода дополнения Шура в программной модели CUDA+OpenMP

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

1 s.kopysov@gmail.ru, 3 Nedozhogin@inbox.ru

ФГБУН «Институт механики Уральского отделения РАН» (ИМ УрО РАН)

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

Аннотация. Эффективное применение метода дополнения Шура на гибридных (CPU/GPU) архитектурах связано с распределением вычислений между центральным процессором и графическими ускорителями. Показано, что формирование матриц дополнения Шура может эффективно выполняться на графическом ускорителе для матриц, состоящих из нескольких тысяч строк и столбцов. Для решения интерфейсной системы предложен параллельный алгоритм метода сопряженных градиентов с явным предобуславливателем, позволяющий достигать существенного ускорения вычислений на нескольких GPU.

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

Ъъомшс QjrAQnQj

Т. 17, № 4 (57). С. 219-229

ВВЕДЕНИЕ

В настоящее время метод дополнения Шура [1, 2] принято рассматривать как гибридный метод решения систем линейных алгебраических уравнений (СЛАУ) [3, 4], сочетающий преимущества как прямых, так и итерационных методов, и учитывающий различные архитектуры параллельных вычислительных систем. На основе дополнения Шура [5] строятся эффективные предобуславливатели для решения СЛАУ различных видов [6]. Раскрыть преимущества данного метода как гибридного позволяют вычисления на CPU/GPU архитектурах. Эффективное применение графических ускорителей, тем более нескольких, в методе дополнения Шура связано с распределением вычислений между CPU и GPU, и соответствующей декомпозицией матриц.

Появление программного обеспечения для вычислений общего назначения на графических устройствах (GPGPU), такого как CUDA, OpenCL и интерфейса прикладного программирования OpenACC (http://www.openacc.org), позволило на операциях линейной алгебры получать ускорение вычислений в десятки и сотни раз по сравнению с центральным процессором.

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

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

Переход к параллельным вычислениям на нескольких графических ускорителях предполагает использование разных технологий параллельного программирования: CUDA - для вычислений на одном GPU, OpenMP - для распараллеливания вычислений между несколькими GPU внутри одного вычислительного модуля.

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

1. МЕТОД ДОПОЛНЕНИЯ ШУРА

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

Пусть область О разбита на пО непересекающихся подобластей (1).

а = а1 и о2 и ...и опп, Ц П nJ = 0,

"а (1)

гв = иац / за.

1=1

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

"а

С(Г, Е) = и о, (V, Е,), '=1

здесь множество вершина графа V - это множество конечных элементов расчетной сетки, множество ребер графа Е - множество смежных конечных элементов, V i с V - множество конечных элементов, образующих подобласть. Далее полагается, что все подграфы О^Ег) связные, в противном случае система уравнений (2) для подобласти О, распадается на несвязанные системы уравнений.

Узлы расчетной сетки образуют множество

V и условно разделяются на внешние Уя - принадлежат границе области и внутренние УI -связанные с узлами подобласти сетки, соответствующей подграфу О^,Е) Из множества внешних узлов выделяются интерфейсные

V с Ув , связанные с узлами из других подобластей.

Для каждой подобласти О, строятся системы уравнений, причем степени свободы, связанные с внутренними и внешними (граничными) узлами, разделяются:

(

А' Ап

А'

АВ1

А' А'

АВВ у

V..! Л

V ив у

( /Л //

V /в У

где индексы I, В относятся к внутренним и граничным степеням свободы соответственно.

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

$вви~в = ~в, (3)

$вв = XI Авв - Ав1АП Ав

/в =Х1 /в Ав1АП /1

-в^п

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

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

А - Т Тт

А11 = ТпТп,

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

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

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

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

А = Т Т1

АII = ТпТп ■

V 6 у

2.

3.

перемен-

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

А' 1в = Ап Аш, {"в ■"2 ) Сформируем матрицы граничной жесткости подобластей

Л

V = А - А А

$вв = Авв Ав1А 1в ■

п, ■ пк

„2 Л

2

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

7в = /в - Ав1 А-1/1, {п1 ■ пв )

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

I

SBBuB = ~, M-'SBB )< C^l + log^fjj-

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

UI = AII fI ~ A IB UB, (ni )

Здесь щ = m ■ VI , nB = m ■ V1

число внут-

ренних и граничных степеней свободы соответственно, m - число степеней свободы в узле сетки, h - шаг сетки, H - размер подобласти, M -предобуславливатель для дополнения Шура. На пятом шаге Алгоритма 1 приведена оценка обусловленности матрицы SBB, а не вычислительная сложность, которая зависит от выбора метода решения системы уравнения. В данной работе в качестве метода решения СЛАУ использовался метод сопряженных градиентов с различными предобуславливателями. Матрицы SBB, Ajj положительно определены и симметричны.

2. РЕСУРСОЕМКОСТЬ МЕТОДА ДОПОЛНЕНИЯ ШУРА

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

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

Предложенный формат DCSR является разновидностью распространенного формата хранения разреженных матриц - формат CSR (Compressed Sparse Row Storage Format). Для матрицы A в формате CSR выделяются три од-

номерных массива, в которых хранятся ненулевые значения [а- | а- Ф 0, 0 < i, у < Ы}, их столбцовые индексы [у | а- Ф 0, 0 < i, у < Ы} и позиции элементов а0, 0 < i < N в двух первых массивах.

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

Ыпг* = тах г=1

Необходимый объем памяти для ленточного формата - (8(2Ыр + 1)Ы), для формата БС8Я -

^ (12Ж2, + 4),

для формата С8Я - (12Ыт + 4(Ы + 1)) байт. В этом случае полагается, что при хранении вещественной величины используется 8 байт памяти, и 4 байта - для целого числа. В симметричном случае величины (2Ыр +1) в оценках сложности алгоритма (и приведенной далее частоты доступа) заменяются на (Ыр + 1), а Ыт и Ыт^ относятся к элементам треугольной матрицы.

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

Преимущества ленточного формата по сложности доступа и добавления элемента компенсируются необходимым объемом памяти (8(2Ыр + 1)Ы) и частотой доступа к элементам матрицы 0((2Ыр + 1)Ы), вместо 0(Ы«г) для форматов БС8Я и С8Я. Здесь

Np = maxj max{/ i i | aф q}

i=1

N

- так называемая полуширина ленты, зависящая от нумерации неизвестных и уравнений, Ы -размерность матрицы, Ыт - число ненулевых

i=1

элементов в матрице, Nnzt - число ненулевых элементов в i-й строке матрицы.

Как показали эксперименты [9], схема хранения оказывает существенное влияние на время вычислений. Отметим, что для ленточных матриц время формирования дополнения Шура составляет порядка 80 %. Переход на формат DCSR для матриц Abb, Ajb, Ajj в одной и той же задаче дал сокращение затрат в четыре раза. Таким образом, в последовательном варианте наиболее эффективным с точки зрения ресурсоем-кости и алгоритмической сложности представляется использовать метод дополнения Шура с разложением Холесского матрицы Ajj и хранением матриц в сжатом формате.

Обработка строк матрицы, хранящейся в формате DCSR, на графическом ускорителе (GPU) потребует для каждой строки: выделения памяти на GPU, копирования строки и далее, в ядре CUDA, выполнения вычислений над строкой и возвращения результата в оперативную память или кэш CPU. Таким образом, потребуется 2N выделений памяти и копирований массивов размера Nnzj, вместо выделения памяти и копирования двух массивов размера Nnz, поэтому для вычислений на GPU матрица переводится в CSR формат.

3. ПАРАЛЛЕЛЬНЫЙ МЕТОД ДОПОЛНЕНИЯ ШУРА И ЕГО ЭФФЕКТИВНОСТЬ

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

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

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

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

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

Сбалансированное распределение вычислительной нагрузки при выполнении шагов 1-4

зависит не от общего числа узлов 1 в подоб-

ластях Qj, а от числа внутренних

V»

Vr

и внешних

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

К

V,

V * ]

приводит к подобластям :

Vr

= 0.

Наглядным примером является задача моделирования напряженно-деформированного со-

стояния пружины, рассматриваемая в данной работе.

Геометрия пружины аппроксимируется неструктурированной расчетной сеткой (рис. 1) с ячейками в виде тетраэдров, число ячеек \У\ =

=174264, число узлов V = 40743. Для получе-

ния подобластей Q

> 0, дуальный граф

полагался взвешенным 0(¥,Е, Ж), с множеством весов Ж = {юк}. Если хотя бы один из узлов конечного элемента не лежит на дО, то вес соответствующей вершины графа юк = 3, в другом

случае о)/, = 1.

Рис. 1. Исходная сетка

б

Рис. 2. Разделенная сетка: а - на 16 подобластей; б - на 1024 подобласти

Проведенные вычислительные эксперименты с различным числом подобластей (пп =16, 32, ..., 1024) показали, что увеличение числа подструктур приводит к увеличению размера матрицы дополнения Шура. Отметим, что число подобластей увеличилось в 64 раза, а размер матрицы 8ВВ только в 1.44 раза N = 66030 в случае пп = 16, и N = 95523 - для пп = 1024). Вместе с тем число ненулевых элементов матрицы дополнения Шура уменьшилось в 18 раз (с Nnz ~ 2,7-108 в случае 16 до Nnz ~ 1,5-107 -для 1024 подобластей), как следствие, уменьшилось ее заполненность с 6.11 % до 0.16 %.

В среднем примерно каждый шестнадцатый элемент в строке SBB является ненулевым (4041 из 66030), при nn = 16 подобластей (рис. 2, а), и примерно каждый шестисотый (154 из 95523) -при разделении на 1024 подобласти (рис. 2, б).

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

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

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

Рассматриваемый в работе алгоритм вычисления обратной матрицы состоит из решений матричной системы вида ЛЦХ = E, где E - единичная матрица n х щ. Система эффективно решается на GPU предобусловленным алгоритмом сопряженных градиентов (см. следующий параграф). Если в правой части системы вместо матрицы E брать Л1В, то ее решением будет матрица

A IB = AiiArn ■

Такое представление позволяет заменить операции обращения матрицы и матричного произведения на решение nB систем, каждая из которых решается независимо, и использовать одновременно несколько GPU. Дополнение Шура SBB или матрица граничных жесткостей вычисляется по соотношениям (3), где ABB е R"b х"в,ЛД е Rnx"I , ABI е RnBxni,A1B e Rn^ ■

Пусть nj - количество столбцов матрицы ЛВВ, пересылаемых на i-й GPU. Каждый графический ускоритель решает nj систем вида:

Aiia = aiB,

здесь a\B - k-й столбец матрицы Л1В,

I е

i . i+1 £ ~B , £ nB,

A'ib = {a1, a 2„

_]=0 ]=0

При реализации независимого решения систем уравнений на нескольких графических ускорителях используется технология ОрепМР. Для этого создаются несколько нитей (число которых равно количеству доступных устройств

i

а

n

в

a

GPU), определяется номер нити и каждой назначается графический ускоритель с тем же номером.

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

(A' 1В У е R"' *, (A' 1В У ={ak\k е

i . i+1 Z , Z ~B ,

j=o j=0

Алгоритм 2 (параллельный алгоритм формирования дополнения Шура на каждой подобласти О/ (индекс 1 опущен): 1. Решаем систему

AIIA IB = AIB ;

{матрицы хранятся на GPU в CSR формате, решение - по столбцам}

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

SBB = A ВВ ~ A BIA IB ;

{SBB и результат произведения матриц -построчно, ABB - CSR формате}

3. Формируем вектор правой части

fB = fB ~ ABIХ;

{x - решение системы AII x = fI}

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

SBBUB = fB;

{SBB формируется в DCSR на CPU, а затем

копируется в CSR на GPU}

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

Uj — x .A. J-Q и^.

{x и A'IB вычислены ранее, и хранятся на CPU}

Это позволяет без дополнительных коммуникаций выполнить оставшиеся операции (произведение и разность матриц, см. соотношение (3)) для каждой матрицы (A'IBy независимо на нескольких GPU, которые используются при вычислении матриц SlBB. Далее формируется глобальная матрица дополнения Шура SBB (шаг 4 в Алгоритме 2), состоящая из локальных SBB ,

принадлежащих i-й подобласти.

Локальные матрицы дополнения Шура хранятся в несжатом формате. Матрица SBB после формирования преобразуется из несжатого к тому формату, в котором будет наиболее удобной с ней работать на GPU, например CSR.

В табл. 1 приведено время параллельного формирования матрицы дополнения Шура в зависимости от числа подобластей и GPU. Во второй колонке представлены данные для последо-

вательного алгоритма, выполняемого на центральном процессоре (Алгоритм 1).

Ускорение параллельного алгоритма (Алгоритм 2), реализованного на GPU, по отношению к CPU, определим как s(np)CPU = tCPU/t(np)GPU, где tCPU - время выполнения последовательного алгоритма, реализованного на CPU, а t(np)GPU -время выполнения параллельной реализации на np GPU. Аналогичным образом введем ускорение s(np)GPU = t(l)GPU/t(np)GPU. В рассмотренных вариантах максимальное ускорение составляет s(8)GPU = 2.7 и достигается при числе подобластей nn = 16, при этом в каждой подобласти находится в среднем по 11000 ячеек сетки.

Формирование матрицы дополнения Шура на двух GPU приводит к ускорению s(2)GPU >1.5, в зависимости от числа подобластей. Использование восьми графических ускорителей дает наименьшее время выполнения для реализации на GPU, но для задач с небольшим числом ячеек в каждой подобласти (< 2500) CPU-реализация оказывается эффективнее. В этом случае при решении большого числа систем уравнений малой размерности для каждой подобласти требуются время на копирование данных между GPU и CPU, которое становится больше, чем время решения систем. В рамках одного GPU (при использовании нескольких) затраты на решение сокращаются, но не покрываются затрат на инициализацию и копирование.

3.2. Решение интерфейсной системы уравнений на GPU

Матрица дополнения Шура

SBB е R > S = S BB = {^j }

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

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

Таблица 1

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

Пп CPU 1 GPU 2 GPU 4 GPU 6 GPU 8 GPU

16 26:23.6 31:52.6 19:25.5 13:o9.9 1o:57.6 o9:49.4

32 o8:oo.6 19:33.9 11:o1.8 o6:41.7 o5:14.o o4:26.5

64 o2:25.1 11:18.8 o6:18.2 o3:44.8 o2:54.7 o2:29.o

128 - - o3:57.2 o2:25.8 o1:58.5 o1:42.5

256 - - o3:14.o o2:o1.o o1:37.7 o1:26.o

512 - - o2:48.3 o1:45.7 o1:26.o o1:15.4

1o24 oo: 18.7 o3:53.4 o2:24.o o1:32.2 o1:17.5 o1:o8.o

В этом случае дополнительные операции, вызванные переходом к предобусловленной системе, сводятся к матрично-векторному произведению zk+i = Mrk+i.

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

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

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

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

1. S,M е RNxN

{M формируется на GPU, матрицы хранятся в CSR формате}

2. u, r,p,q,z е RN

{вектора хранятся в памяти GPU, копии на CPU нет}

3. ro ^ f

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

4. u0 ^ 0

{инициализация выполняется на GPU}

5. zo ^ Mroo {выполняется на GPU}

6. po ^ Zo

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

7. po ^ (ro, zo)

{здесь и далее выполняется

с помощью функции cublasDdot}

8. while ||rI-||2/||b||2 > s do

9. qi ^ Spi

{выполняется на GPU или по формуле

(4)}

Ю. ai ^ (ri, zi)/(qi, pi)

{вычисляется с помощью cublasDdot}

11. ui+i ^ ui + a pi

{операция выполняется с помощью функции cublasDasxpy}

12. ri+1 ^ ri - a qi

{операция выполняется с помощью функции cublasDasxpy }

13. Zi+1 ^ M ri+1 {выполняется на GPU}

14. pi+1 ^ (ri+1, Zi+1)

{вычисляется с помощью cublasDdot}

15. Pi+1 ^ pi+1/pi {вычисляется на CPU}

16. pi+1 ^ Zi+1 + Pi pi {последовательное использование функций cublasDscal и cublasDasxpy}

17. end while

В вычислениях каждой координаты вектора результата используется от 2 до 32 потоков в зависимости от разреженности матрицы. Пре-добуславливатель вычисляется на GPU или на CPU в зависимости от его типа.

3.3. Особенности решения интерфейсной системы на нескольких GPU

Для решения интерфейсной системы (3) реализован блочный алгоритм метода сопряженных градиентов с распределением вычислений между np графическими ускорителями

средствами OpenMP. В этом случае Алгоритм 3 выполнялся в параллельной области, созданной директивой parallel. Результаты скалярных произведений в отдельных нитях OpenMP помещаются в общие переменные и суммируются при помощи директивы atomic, с последующей барьерной синхронизацией.

При разделении на блоки, матрица S = {sj} системы (3) представлялась графом GS(V, E), где V = {i} - множество вершин графа, образованное строчными индексами; E = {(i, j)} - множество ребер графа, образованное парами строчных и столбцовых индексов ненулевых элементов матрицы S; число вершин графа равно размерности S. Граф GS разделяется на np частей многоуровневым алгоритмом [11]. На основе вычисленного разделения, каждой вершине графа ставится в соответствие номер графического ускорителя k. Исходя из назначенного номера GPU, вершины графа подразделяются на внутренние и граничные (связанные хотя бы с одной вершиной, которой сопоставлен другой номер графического ускорителя).

На основе полученного разделения в каждом блоке Sk выделялось несколько матриц: Sk ]

k - матрица, элементы которой связывают

S['k ,bk ] Sb ''k ] вершины в блоке; Sk , Sk - матрицы, связывающие внутренние вершины с граничными;

S[bk 'bm ]

Sk - матрица, связывающая граничные вершины k-го блока с граничными вершинами

m-го блока. Здесь k Ф m и число блоков.

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

k, m е[1, np ]

, где np -

f

S -

S[ 'i' 'i] S ['i'bi] S [bi ''i] s [bi,bi]

0

A1

Si

0

[bl'bn„ ]

o['np np ] Sn„

¡¿X'rip 'bnp ] Sn

0[bnp ''np ] ЛКр 'bnp ] Sn„ Sn„

0 0

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

РТ = (, Р1Р'к, РЬк >•••Р'Пр, РЬПр )> на каждом ОРИ вычисляются два вектора:

nb _ 0[bk ''k ]„' I V V[bk 'bm ]„b .

qk = Sk Pk + Л Sk pm;

m=1

(4)

4k

i - S[!k'!kp + s[k'bk]pb

Pk,

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

Разработанное программное обеспечение для матрично-векторного произведения представлено в Листингах 1-2. Этот программный код является частью параллельной области программы метода сопряженных градиентов и

предназначен для вычисления q' и qb графическим ускорителем k. Барьерная синхронизация обеспечивает получение векторов к моменту начала вычислений на CPU (листинг 1, строка 11) и GPU (там же, строка 25).

При выполнении матрично-векторного произведения, текстурная память привязывается к полю V структур d_p_b и d_p_i в функции MatrVectVul (Листинг 2, строка 24 и 44). Собственно произведения вычисляются kernel-функцией dev matrVectMul, в которой используется текстурная память - для координат вектора p и разделяемая - для доступа к элементам матрицы и сбора промежуточных результатов.

Из приведенных в табл. 2 результатов видно, что решение интерфейсной СЛАУ методом сопряженных градиентов требует меньше временных затрат в алгоритме, использующем GPU. Получено ускорение s(1)cro = 251 для nQ = =64 и s(8)Gro = 3.5 для nQ = 16. С увеличением числа подобластей количество ненулевых элементов в полученной матрице дополнения Шура уменьшается - это приводит к тому, что эффективность использования нескольких GPU для решения интерфейсной системы уменьшается. Например, при разделении на 1024 подобласти s(8)GPU = 1.3.

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

Листинг 1. Фрагмент параллельной области программы метода сопряженных градиентов

1 cudaMemcpyAsync(h.p.b[cpu_thread_id].V, d-p-.b.Vj sizeo f (double) *d_p_b . dim

i cudaMerocpyDeviсeToHost ,streaml); a // вычисление

j MatrVectMul(maxDimGrid , d_Ap_i > d_p_i, d_p_b, d.Ai ; d^Aib, straam2) ;

it // вычисление S.k~{[b-k,i_k]}*p~i-k + Jc [b.k , t>_ к ]} *p " !>„Jc

ft MatrVactMul(roaxDimGrid„ d_Ap_b , d.p.b , d.p.i , d_Ab, d^Abi , 3traam2);

т

a // вычисление g~b_A: \ sum. {m= J } '{irK'n-p } S.k~{[b_k,b.k]}*p~b_k

n for(i-0;iih.Ap.b.dim;it+) h.Ap.b,V[i]"0; lft -¿(pragma о mp barrier

ii fo г ( i=0; i < виbdamain ; i

19 { _

i л if(i) = сpu_thread,id) и {

in for(j=0; j < Abb [i ] . dim ; j+0

lft for(iut 1-Abh[ij . NL[j] ; 1< АЬЪ [1] . NLM +1] ; 1 + +)

ii {

i* h_ Ap_ b . V [j ] + = AbbCi] . V[1]*h_p_b [i] . V [Abb[i] ♦IfC [1]] i

м- }

in }

ii }

11 cudaMemcpy Async (d_ Ap.b.loc . V , h.Ap.b.V, si r.vtt 1' ( climhlr ) * h A b , dim ,

14 cudaMemcpyHoetToDavice , itr(ial)j H ¿iprajjiiia omp barrier

in cublseOaxpy(d.Ap.b,din, 1,0, d_Ap.b.lос ,V, 1, d.Ap.b.V, 1);

Листинг 2. Функция, вычисляющая вектор q'k

i void MatrVectMul (CDtUt tut HiaxDimGr id ,

i Vector d_Ap, Vector d,p.i, Vector d^p_b ,

я CSRmatrix d.Ai j CSRmatrix d_Aibj

.1 cudaStream.t stream)

i {

tt const size.t dimBlock = 123;

т // заполнение вектора q нулями

ч vec_def <<< maxDimGrid, dimBlock >>> (d_Ap.Vj d_Ap,dim);

о

lft iut nnz_per_row - d.AiLMumEl I d.Ai .dim;

ii unsigned iut thr.per.TBc;

li

ii iF(nnz.per.row 2) thr.per_vec=2;

и else iF(nnz.per.rou <= 4) thr_per_vec=4;

in else iF(nnz.per.row <= 3) thr.per.vec-3;

lft else if(nnz.per.ro» <- 16) thr.per_vec=16;

U else thr.per_vec=32;

m aiza.t vaca.per.Ы = dimBlock / thr.psr.vsc;

90

91 aize.t DimGrid = std ; : nin <tut >( maxDimGrid ,

ii (d.Ai,dim +■ iraqa .per _bl - l) / vacs _per _ bl ) \

It

ii cudaBindTaxtura(0, tax.cgb, d.p.i,V, и! иео f (double) # d.p. i , dim) ;

in // о i.i ч а сл с hi a e лераы ас слагаемых из соотношений N

lft dev.MatrVectMul <•; < DiraGrid , dirafll ос k ,0 , str ват > > > (d.Ap.V, d.AKlf,

ii d_Ai,NC, d_Ai.NL, d_Ai,dim,vecs_por_Ы , thr_por^uaс>;

i.i cudaUnbindTextur&(t ex_egb);

:ni nni.pec.ro» " d_Aib>Num£l / d.Aib.dimi

31

ft i f ( nni_ ре r_ r ow <= 2) i.tir_ per ec = 2;

a elst> IF (aiz. ре r_ row <= 4) i.hr.per_vее =4

i f ( ппг. ре r_ r ow <= B) i.tir_ per ec = 9 j

44 else i f ( nns_ ps r_ r ow <= 16) i.tir_ per ec = 16 \

м olitii thr_per_wec = 32 ( ii

it vecs_per~bl = dimSlock / thr_pcr_vsc;

hi DimGrid=Btd::min<iul. >(max DimGrid ,

ii (d.AiЬ.dim + veci_per-1)/veiа_раг_bl);

ii

14 eudaBindTextur-H (0 , tex.cgb, d_p_b,V, ai lieof (double) + d_p_b. din) j

л // вычисление вторь гз сляг ясти uj Со отношен ий (¿)

т dsif.HatrVectHul < < < Di hiGrid , di mfll ae k , 0 , e tr ваш > > > (d_Ap . V , d_ Aib , V ,

hi d_ Aib . HC , d_ Aib . NL , d_ Aib . dim , vbc a_per_bl , thr _per _ vsc )

ir CUdaUnbindTeXturв(tegb);

14 }

Таблица 2

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

Па CPU 1 GPU 2 GPU 4 GPU 6 GPU 8 GPU

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:01: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

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

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

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Haynsworth E. V. On the Shur complement // Basel Mathematical Notes. 1968. № 20. 17 p.

2. Przemieniecki J. S. Theory of Matrix Structural Analusis. New York: McGaw-Hill, 1968. 480 p.

3. Giraud L., HaidarA., Saad Y. Sparse approximations of the Shur complement for parallel algebraic hybrid solvers in 3D // Numerical Mathematics. 2010. Vol. 3. P. 276-294.

4. Rajamanickam S., Boman E. G., Heroux M. A. ShyLU: A hybrid-hybrid solver for multicore platforms // IEEE 26th Int. Parallel and Distributed Processing Symposium (IPDPS), 21-25 May 2012. P. 631-643.

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

6. Корнеев В. Г., Енсен С. Эффективное предобуслав-ливание методом декомпозиции области для р-версии с иерархически базисом // Известия вузов. Математика. 1999. Т. 444, № 5. С. 37-56.

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

8. Kopysov S. P., Krasnopyorov I. V., Novikov A. K., Rychkov V. N. Parallel distributed object-oriented framework for domain decomposition // Domain Decomposition Methods in Science and Engineering. Springer, 2005. Vol. 40. P. 605614.

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

10. Копысов С. П., Новиков А. К., Сагдеева Ю. А. Решение систем уравнений метода Галеркина с разрывными базисными функциями на графическоем ускорителе // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2011. № 3. С. 137-147.

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

ОБ АВТОРАХ

Копысов Сергей Петрович, проф., зав. лаб. выч. и инф. технологий. Дипл. инж.-мех. (Ижевск. мех. ин-т, 1988). Д-р физ.-мат. наук (ИММ РАН, 2007). Иссл. в обл. параллельных вычислений.

Кузьмин Игорь Михайлович, мл. науч. сотр. той же лаб. Дипл. магистр математики (УдГУ, 2009). Иссл. в обл. параллельных вычислений.

Недожогин Никита Сергеевич, асп. той же лаб. Дипл. математик, сист. программист (УдГУ, 2011). Иссл. в обл. параллельных вычислений.

Новиков Александр Константинович, ст. науч. сотр. той же лаб. Дипл. инж.-мех. (ИжГТУ, 1994). Канд. физ.-мат. наук (ИММ РАН, 2005). Иссл. в обл. параллельных вычислений.

METADATA

Title: Parallel algorithms of the Shur complement method in

software model CUDA+OpenMP. Authors: S. P. Kopysov, I. M. Kuzmin, N. S. Nedozhogin, and A. K. Novikov.

Affiliation: Institute of Mechanics Ural Branch of Russian

Academy of Sciences , Russia Email: nedozhogin@inbox.ru. Language: Russian.

Source: Vestnik UGATU (scientific journal of Ufa State Aviation Technical University), vol. 17, no. 4 (57), pp. 219-229, 2013. ISSN 2225-2789 (Online), ISSN 1992-6502 (Print). Abstract: Implementation of the Shur complement method on the hybrid (CPU/GPU) architecture is effective because of quality of computing distribution between CPU and GPUs. We have show that the formation of the Shur complement matrix can be performed efficiently on the GPU for matrix consisting of several thousand rows and columns. To solve the interface system, we proposed parallel algorithm of the conjugate gradients method with explicit preconditioning. This helps to achieve significant acceleration of the computing on several GPUs. Key words: Shur complement method; parallel algorithms;

conjugate gradient method; GPU. References (English Transliteration):

1. E. V. Haynsworth, "On the Shur complement," Basel Mathematical Notes, no. 20, 1968.

2. J. S. Przemieniecki, Theory of Matrix Structural Analusis. New York: McGaw-Hill, 1968.

3. L. Giraud, A. Haidar, and Y. Saad, "Sparse approximations of the Shur complement for parallel algebraic hybrid solvers in 3D," Numerical Mathematics. vol. 3, pp. 276-294, 2010.

4. S. Rajamanickam, E. G. Boman, and M. A. Heroux, "ShyLU: A hybrid-hybrid solver for multicore platforms," IEEE 26th Int. Parallel and Distributed Processing Symp. (IPDPS), 2012, pp. 631-643.

5. D. K. Fadeev and V. N. Fadeeva, Computational Methods of Linear Algebra, (in Russian). Moscow: Fizmatgiz, 1960.

6. V. G. Korneev, S. Ensen, "Efficient pre-domain decomposition method for the p-version with a hierarchical basis", (in Russian), Izvestiya Vuzov. Matematika (Proc. Higher Education. Mathematics), vol. 444, no. 5, pp. 37-56, 1999.

7. S. P. Kopysov, I. V. Krasnopyorov, and V. N. Rychkov, "Object-oriented domain decomposition method", (in Russian), Vychislitel'nye Metody i Programmirovanie (Numerical Methods and Programming), vol. 4, no. 1, pp. 176-193, 2003.

8. S. P. Kopysov, I. V. Krasnopyorov, A. K. Novikov, and V. N. Rychkov, "Parallel distributed object-oriented framework for domain decomposition," Domain Decomposition Methods in Science and Engineering, vol. 40, pp. 605-614, Springer, 2005.

9. S. P. Kopysov, "Optimal separation of the parallel substructure method," (in Russian), in Proc. 5th All-Russian seminar on "Grid methods for solving value problems and applications", Kazan: KGU, 2004, pp. 121-124.

10. S. P. Kopysov, A. K. Novikov, and J. A. Sagdeeva, "Solving simultaneous equations Galerkin method on a discontinuous basis functions on GPU," (in Russian), Vestnik Udmurtskogo Universiteta. Matematika. Mechanika. Kompyuternye Nauki (Bulletin of Udmurt State University. Mathematics. Mechanics. Computer Science), no. 3, pp. 137-147, 2011.

11. G. Karypis, V. Kumar, "Parallel multilevel k-way partitioning scheme for irregular graphs," SIAM Rev., vol. 41, no. 2, pp. 278-300, 1999.

About authors:

Kopysov, Sergey Petrovich, Doctor of Physics and Mathematics, Professor, Institute of Mechanics, Ural Branch of Russian Academy of Sciences.

Kuzmin, Igor Michailovich, Junior Researcher, Institute of Mechanics, Ural Branch of Russian Academy of Sciences.

Nedozhogin, Nikita Sergeyevich, Post-graduate student, Institute of Mechanics, Ural Branch of Russian Academy of Sciences.

Novikov, Alexander Konstantinovich, Candidate of Physics and Mathematics, Senior Researcher, Institute of Mechanics, Ural Branch of Russian Academy of Sciences.