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

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

Том 16, № 4, 2011

Параллельный алгоритм реализации метода встречных циклических прогонок для двумерных сеточных областей*

Л. В. Логанова1, Д. Л. Головашкин2 1 Самарский государственный аэрокосмический университет, Россия 2Институт систем обработки, изображений РАН, Самара, Россия e-mail: tk@smr.ru, dimitriy@smr.ru

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

Ключевые слова: параллельные алгоритмы, метод встречных циклических прогонок.

Введение

Применение параллельных вычислений в математическом моделировании открывает новые возможности для исследований физических явлений, позволяя рассматривать изучаемые процессы на более длительных временных интервалах либо на более обширных областях пространства. В частности, с появлением доступных многопроцессорных систем развитие компьютерной оптики [1] получило новый импульс. Так, современный переход к нанооптике обусловливает популярность Finite-Difference Time-Domain (FDTD) метода моделирования распространения электромагнитного излучения в рамках строгой теории дифракции. Известная монография Yu Wenhua и др. [2] подытоживает многолетнюю практику составления параллельных алгоритмов FDTD-метода. Однако в ней остались не рассмотренными важные разностные схемы Zheng, Chen, Zhang [3] и параллельные алгоритмы для этих схем.

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

* Работа выполнена при финансовой поддержке РФФИ (гранты 10-07-00453-а и 10-07-00553-а).

Кроме того данный алгоритм рекомендуется использовать при изучении других задач вычислительной физики с периодическими краевыми условиями, решение которых приводит к системам линейных алгебраических уравнений (СЛАУ) ленточного вида.

Параллельный алгоритм основан на методе циклических прогонок [4] и является развитием подходов, предложенных в работах [5, 6] для метода встречных прогонок. Известны и другие прямые методы решения подобных СЛАУ на многопроцессорных системах: декомпозиции области и циклической редукции. Однако, как показано в [7], они менее эффективны, чем алгоритм встречных прогонок, либо по числу арифметических операций (декомпозиция области), либо по объему коммуникаций (циклическая редукция),

1. Метод встречных циклических прогонок для двумерных сеточных областей

Положим, что для решения некоторой задачи математической физики (например, теплопроводности или электродинамики) используется двумерная сеточная область wh = х w2h2, где wlhí = {x!(i) = ihui = 0, l,...,M - 1}, = {x2{j) = jh2,j = 0,1,...,N — 1} Ha wh задана сеточная функция gij = g(ihi,jh2). После составления некоторой разностной схемы (либо расщепления, либо переменных направлений [3, 8]) необходимо найти решение совокупности СЛАУ: каждой строке и столбцу сеточной об-wh

задача из [4]

—ао ум-i + соуо — boVi = fo, i = 0, —ai y i-i + CiVi — biVi+i = fi, 1 < i < M — 1, Ум = Vo.

В целом решение всех разностных уравнений схемы (на одном временном слое, если

M N wh N

M

Ak Ук = bk, (1)

где Ак е RMxM, г/к, Ък е RM, ук = (дк,о,-,9к,м-1), к = 0: N -1, aij ф 0 для i < j + 1, j < i + 1, а также aiM = 0 амi = 0,

Aiyi = bt, (2)

где Ai G RNxN, yhbi e RN, y¡ = (g0,i,gN-i,i)T, l = 0 : М - 1, a¿,- ф 0 для i < j + 1, j < i + 1, а также aiN = 0 aNi = 0,

Для решения СЛАУ (1) запишем метод циклических встречных прогонок. Как и в методе циклической прогонки [4], представим решение в виде линейной комбинации сеточных функций

Vi = Ui + yoVi, 0 < i < M — 1. (3)

Для отыскания щ, vi обратимся к методу встречных прогонок [4]. Тогда

а23 0 bi a2i

«2 = -, Р 2 = -, 72 = -, (4)

a22 a22 a22

_ а(г+1)(г+2) п _ Ьг + а(г+1)г А

СИг+1 — -, Рг+1

а(г+1)(г+1) — а(г+1)г аг а(г+1)(г+1) — а(г+1)гаг

«(г+1)г7г

'г • О-Т / [- ч

7г+1 =--—--, г = 2, га — 1, (5)

а(г+1)(г+1) — а(г+1)г аг

. аМ(М-1) ЬМ-1 , аМ 1 /„л

(м-1 = -, Г]М-1 = -, уМ-1 = -, (6)

аММ аММ аММ

. _ «г(г-1) _ Ьг-1 + «г(г+1)Пг

Цг-1 — -X, ^г-1 — -

агг — аг(г+1)Цг агг — аг(г+1)Цг

ф._1 =-; г = М — 1, га (7)

агг — аг(г+1)Сг

«г = «¿+1^+1 + А+Ь ^г = «г+1^+1 + 7г+Ь г = ТО- 1,1, (8)

«г+1 = (г+1П + Пг+1, ^+1 = (г+1^г + "г+1,

г = га, М — 2, г»м = 1, «м = 0, м0 = 0, г>о = 1, (9)

Ьо + Й1М ММ-1 + ^12^1 Уо =-,

а11 — а1М ^М-1 — а12^1

где а,в,7,С,П,"_ прогоночные коэффициенты, у,п — промежуточные сеточные функции, Значения искомой сеточной функции уг найдем по формуле (3), СЛАУ (2) решается аналогично,

2. Параллельный алгоритм для четырехпроцессорной вычислительной системы

Приступая к синтезу четырехпроцеееорного параллельного алгоритма, проведем двумерное разбиение сеточной области и^ на четыре подобласти (рис, 1) и^, и^и^и^ (где припадлежит г-му процессору). На первом этапе алгоритма процессоры 1, 3 выполняют встречные прогонки для строк с 1 по N/2, процессоры 2, 4 — для строк с N/2 + 1 по N, производя вычисления по (4)-(7) и обмениваясь между собой прогоночны-ми коэффициентами (рис, 1, о) (например, процессор 1 производит обмен с процессором

3, процессор 2 — с процессором 4), Обратный ход сопровождается определением промежуточных значений сеточных функций для тех же строк по формулам (8), (9) и обменами этими значениями (рис, 1, б), В завершение процессоры 1 и 2, получив граничные значения промежуточных сеточных функций от процессоров 3, 4 соответственно, вы-

уо

сеточпых функций согласно (3) выполняют все процессоры одновременно (рис, 1, в). На втором этапе процессоры 1, 2 и 3, 4 осуществляют встречные прогонки для столбцов с 1 по М/2 и с М/2 + 1 по М, вычисляя значения прогоночных коэффициентов по (4)-(7) и обмениваясь ими (рис, 1,г), Например, процессор 1 производит обмен с процессором 2, процессор 3 — с процессором 4, Во время выполнения обратного хода определяют значения промежуточных сеточных функций согласно (8), (9) для тех же столбцов (рис, 1,д). После этого процессоры 1, 3, получив граничные значения промежуточных

уо

значение процессорам 2, 4, Выполняя прямой ход, согласно (3) процессоры определяют решение (2) (рис, 1,е),

а

4 4

2 4 -► 4 <-

г

4 у г 4 Г

к 4 к к к 4

б

4 4

г

4 4 4

д

4 к к к 4

2 4 г у 4 Г

в

4 4

4 4

г

е

У Г У 4 г

к 4 к к к 4

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

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

3. Параллельный алгоритм с циклическим разбиением

Распространим идею алгоритма с четырех процессоров на Р. Описание алгоритма будет сопровождаться изложением общего случая (Р процессоров) и примером для 16 процессоров. Выполним циклическое разбиение двумерной области на Р2/4 подобласти (64 для 16 процессоров), где Р = 29, q е N (рис, 2) .Процессор р (0 < р < Р) па г-м шаге будет производить вычисления для подобластей и>рр(г).

4Д) 4(2) 4(3) 4(4) 4(4) 4°(3) 442) 4(1)

4« 4(2) 4(3) 4(4) 4° (4) 443) 4(2) 4(1)

4« 4(2) 4(3) 4(4) 444) 4(3) 4(2) 4°(1)

4( 1) 4(2) 4(3) 4(4) 4(4) 4(3) 4° (2) 441)

<(1) 4(2) 4(3) 4(4) 45(4) 44(3) 43(2) 42(1)

4« 4(2) 4(з) 4(4) 42(4) 45( 3) 4.4(2) 43(1)

4« 4(2) 4(з) 4(4) 43(4) 4.2(3) 4.5(2) 44(1)

4« 4(2) 4(з) 4(4) 44(4) 43(3) 42(2) 45(1)

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

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

Шаги первого этапа алгоритма, связанного с циклическими прогонками по строкам сеточной области, приведены в табл. 1,

На г-м шаге (см, табл. 1) процессоры (с номерами из столбца 2) производят вычисления по алгоритму (столбец 4) над подобластью, указанной в столбце 3, В столбце 5 определены коммуникации, которые осуществляет каждый процессор. Описание коммуникаций выполнено в нотации Дж, Голуба, предложенной в фундаментальной работе [10], Так, запись вида if p = 0 then send (a, в, Y, P/4 — 1) (столбец 5, шаг 1, процессор 0 < p < P/4 — 1) означает, что 0-й процессор посылает значения a,P,Y, определенные

Таблица 1. Шаги первого этапа алгоритма

Шаг (г) Подоб- Вычис-

Процессоры ласть ляемые величины Коммуникации

1 2 3 4 5

1 0 < р < Р/4 - 1 4« a, ß, Y if p = 0 then send (a, ß, 7, P/4 — 1) eise send (a,ß,^,p — 1)

Р/4<р< Р/ 2 - 1 <(1) а, ß, Y if p = P/2 - 1 then send (a, ß, 7, P/4) else send (a, ß, 7, p + 1)

Р/2 < Р < 3 • Р/4 - 1 <(1) if p = P/2 then send ((, £,tp,S- P/4 - 1) eise send ((, p — 1)

3 • P/4 < р < Р - 1 4« \fp = P-l then send ((, £,tp,S- P/4) eise send ((, £,<p,p + 1)

С 2 по 0 < р < Р/4 - 1 4(0 а, ß, y if p = P/4 — 1 then recv (a, ß, 7, 0)

Р/4 - 1 4(0 a, ß, y eise recv (a, ß, 7, p + 1) if p = 0 then send (a, ß, 7, P/4 — 1) else send (a, ß, 7, p — 1)

Р/4 < р < Р/2 - 1 4(0 4(0 a, ß, y a, ß, y if p = P/4 then recv (a, ß, 7, P/2 - 1) eise recv (a,ß,^,p — 1) if p = P/2 - 1 then send (a, ß, 7, P/4) else send (a, ß, 7, p + 1)

Р/2 < р < 3 • Р/4 - 1 4(0 a, ß, y if p = 3 • P/4 - 1 then recv ((, V, Pß)

4(0 eise recv P + 1) if p = P/2 then send ((, ¿, 3 • P/4 - 1) eise send (£, £,<p,p — 1)

3 • Р/4 < р < Р - 1 4(0 4(0 if p = 3 • P/4 then recv ((, £,<p,P — 1) eise recv ((,£,<p,p — 1) \fp = P-l then send ((, £,tp,S- P/4) eise send ((, £,<p,p + 1)

Р/ 4 0 < р < Р/2 - 1 <(Р/ 4) a, ß, y send (a, ß, + P/2) recv ((,£,<p,p + P/2)

Р/2 < р < Р - 1 <(Р/4) send ((,£,<p,p-P/2) recv (a, ß,Y,p — P/2)

1 2 3 4 5

С Р/4 + 1 0 < р < Р/4 - 1 u, V if p = 0 then recv (u, v, P/4 — 1)

по Р/ 2 - 1 U, V recv (u, v,p — 1) if p = P/4 — 1 then send (u, v, 0) else send (u, v,p + 1)

Р/4<р< Р/2 - 1 U, V if p = P/2 — 1 then recv (u, v, P/4) else recv (-u, v,p + 1)

U, V if p = P/4 then send (-u, v, P/2 — 1) else send (-u, v,p — 1)

Р/2 < р < 3 • Р/4 - 1 U, V if p = P/2 then recv (u, v, 3 • P/4 — 1) else recv (u, v,p — 1)

U, V if p = 3 • P/4 - 1 then send (u, v, P/2) else send (u, v,p + 1)

3 • Р/4 < р < Р - 1 U, V if p = P — 1 then send (u, v, 3 • P/4) else send (-u, v,p + 1)

U, V if p = 3 • P/4 then send (-u, v, P — 1) else send (u, v,p — 1)

Р/ 2 0 < р < Р/2 - 1 U, V yo recv (-U, v,p + P/2) send (yo,P + P/2)

Р/2 < р < Р - 1 <{Р/ 2) U, V send (u, v,p — P/2) recv (y0,v,p - P/2)

С Р/2 + 1 по 3 • Р/4 0 < р < Р - 1 У

по формулам (4), (5), процессору Р/4 — 1, остальные процессоры указанного диапазона посылают значения а, в, 7 процессорам с номерами на единицу меньше своих номеров. На втором этапе аналогично выполняются встречные циклические прогонки по столбцам,

Р

соседями. Для нашего примера процессор 0 во время выполнения прогонок по строкам совершает обмены с процессорами 2 и 8, а прогонки по столбцам сопровождаются обменами с процессорами 1 и 4, Таким образом, топология коммуникаций предложенного алгоритма — тор,

В табл. 2 приведены сравнительные характеристики параллельных алгоритмов встречных циклических прогонок с линейным [9] и циклическим (из настоящей работы) разбиением сеточной области. При этом количество арифметических действий по линейному алгоритму составляет С ■ М, длительность выполнения арифметической операции — та, длительность одной коммуникации между процеееорами — тк, Приведен-

Таблица 2. Сравнительные характеристики параллельных алгоритмов

Вид разбиения сеточной области Линейное Циклическое

Простои Отсутствуют

Коммуникации 2 -P-N-Tk ■И) 3 • (TV + М) • rfc + - • Р • rfc - 2 • rfc

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

4. Результаты вычислительных экспериментов

Для подтверждения эффективности предложенного алгоритма проведены вычислительные эксперименты. Так, при решении классической неявной разностной схемы [4] двумерной задачи теплопроводности методом покоординатного расщепления рассматривалась квадратная сеточная область ^ размером начиная с 4000 и заканчивая 120 000 по каждому измерению. Переход от слоя к к слою к +1 по времени осуществлялся следующим образом:

Un,m Цщт _ д 77 _ 1 fk

^ и п,т — 2 Jn,m

- Aiищт = - Lfcm, п = 1,2,N - 1, т = 1,2,М - 1

с циклическими граничными условиями при n = 0, N,

uht} - и

n,m

n,m

Т

= в = 1,2,..., TV-1, m = 1, 2,..., М — 1

с циклическими граничными условиями при т = 0, М. Здесь — сеточная функция в узле (пк 1, т/г2) на /с-ом временном слое, ищт — изменение II на временном интервале т за счет теплопроводности только по направлению х. При этом

Л к — ОЛп + Пп Л к — ОЛ к + Л к

д ттк _ Лп+1,т 2Лп,т + Лп-1,т д Ттк _ Лп,т+1 2Лп,т + Лп,т-1 Л1лп,т - -^-> 2 п,т ~ -

h2

h2

Эксперименты проводились на кластере лаборатории "Математические основы дифракционной оптики и обработки изображений" Самарского государственного аэрокосмического университета им. акад. С.П. Королева. Использовались четыре ЭВМ AMD Sempron(tm) Processor 3000+ (1.8 ГГц), ОЗУ 1 Гб, связанные локальной сетью Ethernet 100 Мб. На компьютерах установлена операционная система Fedora. Алгоритмы реализованы на языке Fortran с использованием библиотеки MPI.

4

3.5 3

2.5 2

1.5 1

0.5 О

N

4000 6000 8000 10000 12000 20000 40000 60000 120000

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

Результаты исследований позволяют говорить об эффективности предложенного алгоритма, Из графика на рис, 3 видна тенденция к увеличению величины ускорения, достигающего максимальной величины 3,8, Проведенное сравнение с пакетом Pardiso библиотеки Intel®MKL [11] показало значительное преимущество разработанного подхода, которое растет с увеличением числа узлов за счет учета ленточной структуры матрицы.

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

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

[1] СоЙФЕР В. А. Дифракционная компьютерная оптика. М.: Физматлит, 2007. 736 с.

[2] Parallel Finite-Difference Time-Domain Method / Yu Wenhua, Mittra Raj, Su Tao, Liu Yongjun, Y. Xiaoling. Artech House Publ., 2006. 274 p.

[3] Taflove A., Hagness S. Computational Electrodynamics: The Finite-Difference TimeDomain Method: 3nd. ed. Artech House Publ., 2005. 832 p.

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

[5] Миренков H.H. Параллельные алгоритмы для решения задач на однородных вычислительных системах // Вычисл. системы. 1973. № 57. С. 3-32.

[6] Бирюкова Л.Н., Четверушкин Б.Н. О возможности реализации квазигидродинамической модели полупроводниковой плазмы на многопроцессорных вычислительных системах // Матем. моделирование. 1991. Т. 3, № 6. С. 61-71.

[7] Головашкин Д.Л. Параллельные алгоритмы решения сеточных уравнений трехдиагонального вида, основанного на методе встречных прогонок // Там же. 2005. Т. 17, № 11. С.118-128.

[8] Самарский A.A., Вавищевич П.Н. Вычислительная теплопередача. М.: Едиториал УРСС, 2003. 784 с.

[9] Елизарова Т.Г., Четверушкин Б.Н. Применение многопроцессорных транспьютерных систем для решения задач математической физики // Матем. моделирование. 1992. Т. 4, № 11. С. 75-100.

[10] Голув Дж., Ван Лоун Ч. Матричные вычисления М.: Мир, 1999. 548 с.

[11] Компилятор Intel® Visual Fortran Compiler Professional Edition для ОС Windows [электронный ресурс] // Центр ознакомительного ПО Intel, http://software.intel.com/ru-ru/articles/intel-software-evaluation-center/ (дата обращения 14.04.10).

Поступила в редакцию 26 мая 2009 г., с доработки — 11 февраля 2011 г.