Многопроцессорная реализация одного итерационного алгоритма для систем с распределённой памятью Текст научной статьи по специальности «Математика»

УДК 519.615.5

С. А. Ишанов, С. В. Клевцур,

К. С. Латышев, Д. И. Пялов

МНОГОПРОЦЕССОРНАЯ РЕАЛИЗАЦИЯ ОДНОГО ИТЕРАЦИОННОГО АЛГОРИТМА ДЛЯ СИСТЕМ С РАСПРЕДЕЛЁННОЙ ПАМЯТЬЮ

Рассмотрены многопроцессорные реализации а-Р-итерационного алгоритма решения систем пятиточечных разностных уравнений.

Some multiprocessors realizations of the a-fi-iteration algorithm for solution of five-dots difference equations systems are considered.

Ключевые слова: система линейных уравнений, разностная схема, а — Р-ите-рация, пятиточечная схема, мультипроцессорная система.

Keywords: system of linear equations, difference scheme, a — Р-iteration, five-dots scheme, multiprocessors system.

Введение

Использование численных методов в задачах математического моделирования приводит к необходимости решения систем линейных алгебраических уравнений (СЛАУ) больших порядков, что влечет за собой большие затраты вычислительных ресурсов. В работе [1] для решения двухмерных разностных уравнений предложен а — Р-итерацион-ный метод, не требуюшдй априорной информации о границах спектра разностного оператора, малочувствительный к сильным изменениям коэффициентов разностной схемы и обладающий высокой скоростью сходимости по сравнению с остальными методами [2]. Несмотря на указанные достоинства метода, время вычислений заметно возрастает при увеличении числа узлов разностной сетки. Так, например, расчет модельной двухмерной задачи теплопроводности в области размером 500x500 узлов с 20 шагами по времени занимает около 100 с, а решение той же задачи на сетке 2000x500 c 10 шагами по времени на том же компьютере занимает уже порядка 600 с. В настоящее время для уменьшения времени расчетов широко используются технологии параллельных вычислений на многопроцессорных системах.

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

Вестник Российского государственного университета им. И. Канта. 2009. Вып. 10. С. 64—73.

метода была опубликована в работе [3], где приводится оценка ускорения вычислений в реализации для 32-транспьютерной системы.

Коротко рассмотрим а — Р-итерационный метод [1] на примере двумерного уравнения теплопроводности [4]

— + div(Dgrad и) = Q = I (1)

&

с граничными и начальными условиями и|хеС = ф(х, Ь), и|( =0 = и0(х).

В прямоугольной системе координат уравнение (1) в случае стационарного коэффициента теплопроводности принимает вид

/ 2 2 йи(х, у, í) _ п| й и(х, у, і) + й и(х, у, í)

йх йу

•_ В

йі

+ ¥(х, у, і),

и1х=0 = Ф0(*), и1х=й =ф« V, и\у=0 =У0«, и|у=Ь = УЬ (0, Я=0 = и0(х).

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

т+1 --Т ( ит+1 2ит+1 + ит+1 ит+1 2ит+1 + ит+1 ^

иі+1, ] 2иі, ] + иі-1, ] + иі, ]+1 2иі, ] + иі, ]-1

иі,І - иі,і. _ В

Ді

Дх2 Дх2

+ ¥т+1

і,І

Вычисление значений и на временном слое т + 1 при известных значениях и и ¥ на слое т сводится к разрешению пятиточечной схемы, представленной в общем виде уравнением

Ві,пиі,П-1 + Кі,пиі-1,П — Сі,пиі,П + Еі,пиі + 1,П + ^і,«иі,« + 1 + ¥і,п _ 0 (2)

с граничными условиями

и1,п _ф2,п^2,п + У2,п , ,п _фМ-1,пим, - 1,п + УN. -1,п ,

иі,1 _фі,2иі,2 + Уі,2, иі,М _фі,Мп-1иі,N.-1 + УіМ„-1.

Пусть решение исходной системы (2) удовлетворяет соотношениям

иі ,п _ аі + 1, пиі + 1,п +Рі + 1,п , иі ,п _ у і-1,пиі-1,п + йі-1,п , иі,п _ аі,п + 1иі,п +1 + Рі,п + ^ и,п _ Уі,п-1иі,п-1 + йі,п-1,

где а, у, ~, ~, р, й, Р, й — неизвестные коэффициенты прогонки. Метод а — Р-итераций, как уже отмечалось, подробно изложен в работе [1] и мы не будем приводить расчетные формулы подробно. С модификациями метода можно ознакомиться в работах [5; 6].

Метод распараллеливания

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

Рис. 1. Разбиение без перекрытия (Цвн — значение на внутренней границе)

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

В случае перекрывающихся областей расчет начинается с постановки искусственных граничных условий на границах подобластей ab и сд (рис. 2). С этими граничными условиями производится расчет сеточной функции в подобластях I и II. После окончания расчета в области I определяются новые значения Ц на линии аЬ, в свою очередь, из решения задачи в области II определяются значения Ц на линии сд. После этого расчет в области I повторяется с новыми граничными значениями Ц на линии сд, а в области II — на линии аЬ. Процесс повторяется до сходимости итераций между подобластями.

а с

Рис. 2. Разбиение с перекрытием

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

Вариант разбиения без перекрытия свободен от указанных недостатков. Однако в этом случае остается непонятным механизм формиро-

вания граничного условия на внутренней границе аЬ в процессе итераций между подобластями.

Итерационный а—Р-метод дает простую возможность формирования условий на внутренней границе в процессе итераций между подобластями. В самом деле, в ходе итераций внутри каждой подобласти вычисляются не значения сеточной функции и, а значения прогоноч-ных коэффициентов. После окончания расчета в области I на внутренней границе аЬ получаются новые значения коэффициентов а и р. А после окончания расчета в области II на этой границе определяются значения коэффициентов у и й. Полученные значения коэффициентов а и р послужат новыми граничными условиями для области II, а значения коэффициентов у и й — граничными условиями для области I. Построенный таким образом алгоритм нахождения внутренних граничных условий служит основой параллельного варианта.

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

I. Выделение необходимой памяти.

II. Задание начальных условий.

III. Процесс по к

1) задание коэффициентов схемы и граничных условий;

2) а-итерационный процесс:

а) первая половина а-процесса (вычисление а, у, а),

б) вторая половина а-процесса (вычисление а, у, ~),

в) проверка сходимости;

3) Р-итерационный процесс:

а) первая половина Р-процесса (вычисление р, й, в),

б) вторая половина Р-процесса (вычисление р, й, й),

в) проверка сходимости;

4) вычисление значения сеточной функции в момент времени к,

5) сохранение результатов.

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

Рассмотрим более детально работу системы на стадиях Ш.2а, Ш.2б, Ш.3а и Ш.3б. Организация вычислительного процесса на каждом из четырех этапов одинакова, поэтому рассмотрим только первую половину а-процесса (вычисление а, у, а). Для произведения расчетов одному процессору необходимо наличие краевых условий для а или у либо уже

вычисленных значений а и у. Тогда алгоритм работы каждого процессора можно представить в виде схемы (рис. 3).

Рис. 3. Схема работы первой половины а-процесса

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

Сделаем оценку для ускорения и эффективности параллельной реализации первой половины а-процесса. Рассмотрим вычислительный процесс в системе с четным количеством процессоров одинаковой производительности. Отметим, что четность количества процессоров не сильно сказывается на общей оценке, а лишь упрощает дальнейшие рассуждения. Как было сказано выше, для распараллеливания исполь-

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

а = — - 2Ь, (3)

п

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

п

прогоночного коэффициента а, второй — у — потребует —(а + 2Ь) + Ь

единиц времени. Тогда до момента, когда первый процессор закончит вычисление у, а последний — вычисления а (остальные процессоры в

это время будут заняты вычислением а), пройдет еще П (а + 2Ь) + Ь

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

п(а + 2Ь) + а + 5Ь (4)

единиц времени. Подставляя значение (3) в выражение (4), получим формулу для относительного времени, необходимого для расчета прогоночных коэффициентов на одном слое по у:

1

Тп = N(1 + -) + 3Ь. п

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

1

Тп = N(1 + -). п

Для нахождения эталонного времени вычисления прогоночных коэффициентов на одном процессоре полученная формула непримени-

ма, однако вывод ее очевиден, время расчета складывается из времени определения коэффициентов а, у, а, т. е.

Т1 = 3К

Тогда теоретическое ускорение параллельных вычислений для про-гоночных коэффициентов составит

К = 3 -+-, (5)

п +1

а теоретическая эффективность

3

Рп = . (6)

п+1

Формулы (5) и (6) имеют смысл при и ^ 2 и наибо лее точны для четных значений п. Из формулы (5) видно, что теоретическое ускорение для а- и Р-процессов не превосходит 3 и стремится к 3 при увеличении п. Максимальная теоретическая эффективность достигается при п = 2 и уменьшается при увеличении п. К сожалению, а- и Р-процессы составляют основу алгоритма и добиться более высокого ускорения тяжело. Отметим, что достаточно ресурсоемкие этапы III.4 и III.5 позволяют увеличить это значение. Общая же оценка ускорения параллельной реализации сопряжена с большими трудностями. Это связано с неоднородностью вычислительного процесса, однако суммарная оценка эффективности алгоритма может быть получена экспериментальным путем. Для этого рассмотрим следующее уравнение теплопроводности:

м( = Дм + /, и\ =0 = 0,

Г=_______4________________3______________

1 ^ 2 1 , п СЧ2 , п СЧ2

В качестве граничных условий взято м|^ )еС = 0, что соответствует

бесконечно мощному холодильнику на границе области. Задача была решена (рис. 4) на сетке разного размера и с разным числом шагов по времени: 1) 500x500, 20 шагов по к 2) 2000x500, 10 шагов по t.

Рис. 4. Пример распределения тепла (Н — горячая область, С — холодная)

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

120 100 80 60 40 20

0 ....

1 PU 2 PU 3 PU 4 PU

Рис. 5. Результаты задачи размером 500x500, 20 шагов по t

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

800 600 400 200

0 , , , .

1 PU 2 PU 3 PU 4 PU

Рис. 6. Результаты задачи размером 2000x500, 10 шагов по t

На сводной диаграмме (рис. 7) наглядно изображено изменение среднего времени вычисления при увеличении числа процессоров.

□ Тест 1

□ Тест 2

□ Тест 3

500

400

300

200

100

0

Я 1

х “^“500x500 (20) 2000x500 (10)

1 Ри 2 Ри 3 Ри 4 Ри

Рис. 7. Среднее время вычислений.

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

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

500x500 (20) 2000x500 (10)

Рис. 8. Ускорение вычислений

500x500 (20) 2000x500 (10)

Рис. 9. Эффективность вычислений

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

Заключение

Задачи математического моделирования требуют больших вычислительных затрат. Эффективно увеличить производительность вычислительной системы можно только путем увеличения числа вычислительных единиц. Однако это требует модификации исходного алгоритма с учетом новой специфики системы. Для рекурсивных алгоритмов реализация параллельного алгоритма затруднена, в особенности для систем с распределенной памятью. Но даже в этом случае использование многопроцессорной системы может уменьшить время решения задачи в разы, в особенности для больших задач. В данной работе был реализован многопроцессорный вариант для a-ß-итерационного метода, основу которого составляют рекуррентные вычисления, проведена теоретическая оценка для ускорения и эффективности параллельной реализации. Практические результаты исследования показывают, что даже в случае сложного итерационного алгоритма можно получить заметный выигрыш в скорости вычислений за счет использования многопроцессорных систем.

Работа выполнена при поддержке гранта РФФИ № 09-01-00628-А.

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

1. Четверушкин Б. Н. Математическое моделирование задач динамики излучающего газа. М., 1985.

2. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М., 1978.

3. Четверушкин Б. Н. Кинетически-согласованные схемы в газовой динамике: новая модель вязкого газа, алгоритмы, параллельная реализация, приложения. М., 1999.

4. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М., 1999.

5. Клевцур С. В., Латышев К. С., Четверушкин Б. Н. Циклический вариант a-ß-итерационного алгоритма // Дифференциальные уравнения. 1988. Т. 24. № 7. С. 1213-1218.

6. Ишанов С. А., Клевцур С. В., Латышев К. С. Алгоритм a-ß-итераций в задачах моделирования ионосферной плазмы / / Матем. моделирование. 2009. Т. 21. № 1. С. 33-45.

Об авторах

С. А. Ишанов — канд. физ.-мат. наук, проф., РГУ им. И. Канта, e-mail: math@dekan. albertina. ru.

C. В. Клевцур — канд. физ.-мат. наук, доц., РГУ им. И. Канта.

К. С. Латышев — д-р физ.-мат. наук, проф., РГУ им. И. Канта.

Д. И. Пялов — асп., РГУ им. И. Канта.

Authors

Dr S. A. Ishanov — professor, IKSUR, e-mail: math@dekan. albertina. ru.

Dr S. V. Klevtsur — assistant professor, IKSUR.

Professor K. S. Latyshev — IKSUR.

D. I. Pylov — PhD student, IKSUR.