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

Вычислительные технологии

Том 13, Специальный выпуск 4, 2008

Параллельный прямой метод решения разделяющихся краевых задач

Д.В. Кныш

Новосибирский государственный университет, Россия e-mail: [email protected]

A parallel direct method for solution of a separable boundary value problem is considered. The method is on the Fourier transformation (applied twice) along with the method of reduction. Efficiency of parallelization of the proposed method is analyzed for of the case of 3D problems. The programming implementation for distributed and shared memory systems along with series of numerical measurements are performed.

Введение

Алгоритмы решения систем линейных алгебраических уравнений (СЛАУ) высокого порядка с разреженными матрицами специальной блочной структуры, возникающих при аппроксимации многомерных краевых задач с разделяющимися переменными, изучались многими авторами (см., например, [1-3] и приведенную там библиографию). Такие методы, обладающие рекордным быстродействием, имеют самостоятельное значение для многих приложений, а также могут быть использованы для построения эффективных предобуеловленных итерационных процессов, применимых к решению более сложных уравнений математической физики. К числу наиболее экономичных алгоритмов относятся методы быстрого преобразования Фурье (БПФ), циклической редукции, а также итерационные неявные методы переменных направлений с оптимальной последовательностью параметров, использующие особенности спектральных и структурных свойств матриц. В работе [4] рассмотрено применение комбинированных прямых и итерационных методов для решения диекретизированных трехмерных краевых задач с полуразделяющимися переменными.

Целью данной работы является анализ эффективности распараллеливания прямого алгоритма для решения трехмерных смешанных краевых задач. В качестве характерного примера рассматривается задача Дирихле для уравнения Пуассона в кубе:

д2и д2и д2и

+ ^+ ^2=/(*>*>>*)> у,г) € П = [0,1] х [0,1] х [0,1], и\Г = д. (1)

После аппроксимации задачи (1) на кубической сетке с шагом к — 1/^ + 1) получаемая СЛАУ может быть записана в форме

(Аи)^к = [(Ах + Ау + А* )и]ч,к —

— 6ui,j,k иг+1 ,],к иг,^'+1 ,к 1 — fi,j,k, (2)

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

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

ющий дискретизации вдоль одной из осей, а правая часть определяется следующим соотношением:

f f(xi,yj,zk)h2,i,j,k G 2,N-l,

fi,j,k = S 2 , , , (3)

l f (xi,yj ,zk)h2 + g(xi' ,yj' ,zk'),i,j,k G{1,N},i,j,k g{0,n +1}.

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

1. Описание численных алгоритмов

Записывая систему (2) в векторном виде

Аи = / = {/*,,■*}, и = {и^}, А = Ах + Ау + А, (4)

отметим, что в данном случае А,и, / - симметричная матрица и векторы порядка N А

3

Vp,q,r = {j = 4>?jЧ>?; P,q,r = 1,..., N},

(5)

х,р / 2 • Ф?г / 2 . уди / 2 . Ахтг (Л =\ - БШ-, = 4 /- БШ-, =\ - 8111-,

у м + 1 м +1' ^ V ^ +1 ^ +1 V ^ +1 л^ + 1

а соответствующие собственные числа равны Лр'9'г = Л^ + Лу + ЛГ , где

^ = 48ш22(^тту л? = 4йп22(]1ТТ)- А' = 4йп22(лГГ7)- (в)

При столбцовой упорядоченности узлов и соответствующих компонентов векторов и = {иг}, / = {/¿}, I = к+ (г —+ (. — 1)Ж2, к, г,. = 1,..., N матрица А2 имеет блочно-диагональный вид А2 = сИа§{А2:}, где каждый блок А2 порядка N имеет собственные векторы = {(Д'г} и соответствующие собствеппые числа ЛГ го (6). У матрицы А2 — те же собственные числа, но каждое из них имеет кратность N2, а ее собственные векторы, соответствующие совокупности N неповторяющихся чисел Л^, записываются в блочном виде (З = {(з, ¿,.7, = 1,...,N}, Здесь каждый подвектор (З = {(Г} имеет порядок N а индексы г,. указывают номер блока, где расположены ненулевые элементы (которые фактически зависят только от индекса к, то не от г,.). Количество разных собственных векторов ('3 и порядок каждого из них равны N3, Раскладывая векторы / и и то базису з"? мы можем записать

N N

/ = { X) Дз }, и = { X) }, (7)

Г= 1 Г=1

гДе Л\' и ь?7- для каждого г,] представляют собой коэффициенты разложений Фурье

1г, ] "г, ]

"одномерных столбцовых" векторов порядка N по соответствующему базису у*

N N

fi¿ = Ylf,j,r' vij = Ui,j,' r- (8)

k=1 k=1

Формулы (7) и (8) представляют собой реализацию прямого и обратного разложений Фурье соответственно. Вычисление всех коэффициентов Фурье с помощью приведенных выражений требует выполнения арифметических операций. Однако если здесь применить БПФ [1], то объем вычислений сокращается до 2N31о§2^

Схема предлагаемого алгоритма решения поставленной трехмерной задачи подразумевает использование двукратного (по двум направлениям) БПФ и алгоритма редукции дня получаемых вспомогательных трехдиагональных СЛАУ,

С этой целью "двумерные" векторы ьТ ={у?.} и /Т ={для каждого г разложим в ряды Фурье по второй переменной у:

N N

ьг = {£ г $ Л, и*' ? = £ ^ $' *,

' J (9)

/г = { еГ=1 /г^Г}, /г = ЕГ=1 /Г.^.

В результате для векторов и*'Т = {и*'г} Е Кп получаем уравнения

+ ЛуI + Л*гI)ь*'г = /= {/*'?}, я,г = 1,..., N. (10)

Решение каждой из этих СЛАУ N^0 порядка проводим с помощью хорошо распараллеливаемого метода редукции, который кратко представим в следующей интерпретации. Множество узлов одномерной сетки Q = {I = 1,..., N} разобьем на 2Р + 1 подобласти 0,р, из которых нечетные называются разделителями и состоят из одного узла, а четные области содержат по Ы узлов, при чем 1 + (Ы + 1)Р = N (см, рисунок). Соответствующие уравнения могут быть записаны в виде

(Ли)р = -ЛрЩ-1 + Брир - Срир+1 = /р, р =1, 2,..., Р, Л1 = Ср = 0, (11)

где ир с нечетными номерами есть скаляры, а с четными — векторы Ы-го порядка. Исключая из (11) подвекторы ир±1 с помощью соотношений

ир+1 = Бр+11(/р+1 — Ар+1ир — Ср+1ир+2) = ^р+1 — ^р+1ир — ^+1^+2-, Бр+1Ьр+1 = [1,..., 0]та, Бр+1Ьр+1 = [0,..., 1]тс,

р

арир-2 + Ьрир + Срир+2 — /р, С1р = АрВр\Ар_ 1, Ьр = Вр — АрВр\Ср_1 — СрВр^Ар+1, ср = СрВр^Ср^.

(12)

(13)

i-2 ¡-1 i ¡+1 i+5

Схема декомпозиции области в методе редукции

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

2. Распараллеливание вычислительного процесса

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

а) выполнение прямого быстрого преобразования Фурье по оси z;

б) выполнение прямого БПФ по оси у;

в) решение N2 трехдпагональных СЛАУ методом редукции по оси x;

г) выполнение обратного преобразования Фурье по оси у;

д) реализация обратного БПФ по оси z, результатом чего является искомое решение.

Проведем одномерную декомпозицию области на непересекающиеся подобласти Q = p

= (J Qt с помощью горизонтальных плоскостей x = X0,..,XP, Qt = {Xt-1 < x < t= i

Xt}, t = 1,...,P, образующих равномерную макросеть в том смысле, что в каждой подобласти находится одинаковое число узлов N3/P,

Подвергнем соответствующему разбиению массив правых частей /, а полученные подмассивы отправим на обработку соответствующему процессору р.

После выполнения двукратного преобразования Фурье по направлениям z и у переходим к решению редуцированных подсистем, правые части которых сформированы в ходе выполнения этапов "а" и "б" алгоритма и располагаются в отдельных процессорах. Проанализируем число коммуникационных обменов, требуемое в ходе выполнения этапа "редукция". Необходимость передачи данных соседнему процессору наступает после вычисления w, v и v для подсчета коэффициентов правых частей уравнений, связывающих узлы-разделители, где задейетвуютея околограничные значения w, расположенные в разных подобластях. Этот обмен требует (P — 1) парных операций посылки-получения массива размерности N2, Отметим, что в силу известных соотношений такой обмен работает эффективнее, нежели N2(P — 1)-я операция посылки-получения одного

числа. Следующим шагом вычисляются ikVx, Vy и данные ik, находящиеся во взаимно-

р

очередь коммуникаций между процессорами с целью сборки всех правых частей системы, связывающей узлы-разделители, в некотором едином процессоре, внутри которого она решается. Это действие требует (P — 1) операций посылки-получения N2 чисел. При обратной рассылке вычисленных значений в узлах-разделителях на свой процессор потребуется вдвое больше операций, а именно 2(P — 1) обменов, так как для сборки полного решения трехдиагональных СЛАУ необходимы оба значения в узлах-разделителях, окаймляющих подобласть. Таким образом, всего коммуникационные потери в ходе работы алгоритма составляют

T = 4(P — 1)(то + tcN2), (14)

где тс — время пересылки одного числа, аг0- время задержки одной операции обмена.

Время выполнения всех арифметических операций на Р процессорах можно оценить величиной

^8jV4log2jV , 12 N*\

J-a — Ta I p + p I .

где первое слагаемое в скобках относится ко всем БПФ, второе — к методу редукции решения п2 трехдиагональной СЛАУ, а та есть среднее время выполнения одной арифметической операции.

Реализация алгоритма на архитектурах с общей памятью (МВС-ОП) отличается от описанного способа распараллеливания на МВС-РП в первую очередь оригинальным применением композиции методов решения трехдиагональных СЛАУ, где решение каждой редуцированной подсистемы этапа 3 алгоритма осуществляется с помощью алгоритма встречных прогонок (см., например [2]), в основе которого лежит идея использования для подобластей разных вариантов методов прогонки, осуществляющих вычисления в противоположных направлениях и встречающихся лишь однажды посредине отрезка после прямого хода. Таким образом задача расщепляется на две части, обработка которых может проводиться параллельно отдельными процессорами. Вообще говоря, использование общей памяти здесь не принципиально, обозначенное разделение может быть применено и на МВС-РП, но привлекательно за счет отсутствия коммуникационных потерь и перспективы применения в дальнейшем этого подхода для реализации на гибридных архитектурах, где каждая редуцированная подсистема располагается и обрабатывается на каком-нибудь узле своим МР1-процеееом, а внутри узла части встречной прогонки обрабатываются каждая своим потоком.

3. Примеры численных экспериментов

Проиллюстрируем эффективность распараллеливания описанного прямого метода на МВС-РП и МВС-ОП в применении к решению симметричных систем семиточечных уравнений, аппроксимирующих задачу Дирихле для уравнения Пуассона в единичном кубе на различных кубических сетках. Функции f и g в (1) выбирались из условия, что точное решение u(x,y,z) равно единице. Реализация алгоритмов осуществлялась

Таблица 1. Временные характеристики прямого метода в MPI (двукратное БПФ + метод редукции)

N \ р 1 2 4 8 16

129 0.17 8.19 • 10"2 3.9 • 10"2 1.8 • 10"2 1.8 • 10"2

0.53 0.35 0.28 0.18 9.76 • 10"2

0.16 8.19 • 10"2 4.29 • 10"2 1.8 • 10"2 1.8 • 10"2

ttotal 0.86 0.52 0.36 0.22 0.16

257 1.56 0.77 0.38 0.15 0.12

5.50 2.83 1.40 0.67 0.37

1.5 0.75 0.36 0.18 8.58 • 10"2

ttotal 8.58 4.38 2.17 1.06 0.58

513 - 11.65 5.55 2.86 1.46

- 28.09 12.94 6.14 2.91

- 11.44 5.51 2.87 1.44

ttotal - 51.20 24.03 11.88 5.83

Таблица 2. Временные характеристики прямого метода в ОрепМР (двукратное БПФ + метод редукции * метод встречной прогонки)

N \th 0 1 2

129 0.17 0.17 0.1

0.35 0.39 0.38

0.17 0.17 0.1

ttotal 0.68 0.72 0.59

257 1.56 1.56 0.83

3.32 3.35 2.98

1.51 1.50 0.80

ttotal 6.41 6.45 4.64

513 22.89 22.84 11.84

33.21 31.10 25.28

22.49 22.47 11.64

ttotal 78.61 76.43 48.78

на Фортране, а все вычисления проводились со стандартной двойной точностью. Эксперименты проводились для трех разных сеток со значениями N = 129, 257, 513, В табл. 1 и 2 приведены результаты измерений времени (в секундах) исполнения параллельных версий алгоритма в системе MPI для 1, 2, 4, 8 и 16 процессоров (на кластере с двух процессорными узлами Itanium 2) и ОрепМР для двух потоков (на четырех-процессорном Itanium 2) ради анализа эффективности использования параллельной версии метода встречных прогонок в последовательной версии рассматриваемого метода, Каждая клетка обеих таблиц содержит три значения, соответствующих времени исполнения соответствующего этапа алгоритма, сверху вниз — общее время ППФ(^) и ППФ(?/), редукция, общее время ОПФ(?/) и ОПФ(^) и полное время решения задачи ttotal.

Заключение

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

Список литературы

[1] Самарский A.A., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978.

[2] Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. Новосибирск: НВМ СО РАН, 2000.

[3] Вшивков В.А., Снытников В.Н., Снытников Н.В. Моделирование трехмерной динамики вещества в гравитационном поле на многопроцессорных ЭВМ // Вычисл. технологии. 2006. Т. 11, № 2. С. 15-27.

[4] Ильин В.П., Кныш Д.В. Параллельные алгоритмы решения разделяющихся краевых задач. СПб.: Изд. СПбПГУ, 2008. С. 107-118.

[5] Богачев К.Ю. Основы параллельного программирования. М.: БИНОМ, Лаб. знаний, 2003.

Поступила в редакцию 20 февраля 2007 г., в переработанном виде — 13 мая 2008 г.