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

Раздел III. Алгоритмическое и математическое обеспечение вычислительных систем

Л.П. Фельдман, И.А. Назарова

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

Моделирование реальных физических, технических и других процессов, описываемых системами обыкновенных дифференциальных уравнений (СОДУ), представляет собой обширный класс задач, для решения которых применение высокопроизводительной вычислительной техники не только оправдано, но и необходимо. Об этом свидетельствует знаменитый список проблем “большой вызов”, в котором такие задачи занимают одно из ведущих мест [1]. Как показала практика, наиболее эксплуатируемым способом создания параллельных методов является распараллеливание хорошо исследованных и многократно апробированных последовательных численных алгоритмов [2]. Современный качественный численный алгоритм решения СОДУ обязательно должен содержать механизм управления шагом интегрирования на основе информации о погрешности решения на каждом шаге интегрирования.

В данной статье рассматривается эффективность альтернативных методов оценки апостериорной локальной погрешности на базе параллельных алгоритмов решения СОДУ для многопроцессорных ВС 81МБ-архитектуры с распределенной памятью.

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

Наиболее простой способ достижения поставленной цели - это удвоение или дублирование шага по правилу Рунге. На основе одной и той же формулы Рунге-Кутта порядка р вычисляют два приближения к решению в одной и той же точке

сначала с удвоенным шагом - 2к, а затем независимо, дважды с половинным шагом - к. Для точного решения при переходе из точки хп к точке хп + 2к, у(хп + 2к) получают две аппроксимации у1(1 шаг 2к) и у2 (2 шага к).

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

у(х + 2к) = у1 + (2к)р (р + 0(кр);

у(х + 2к) = у 2 + 2(кр )ф + 0(кр);

Ду = у 2 - у1;

у(х + 2к) = у з = у 2 +-Д— + 0(кр+1).

2 — 1

Процесс уточнения решения с половинным шагом у2, как видно из (1), повышает точность решения на единицу и получил название локальной экстраполя-

ции. Данный способ оценки погрешности достаточно прост, однако имеет большие накладные расходы: объем вычислений на один узел сетки возрастает почти втрое.

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

5

Уп+1 = Уп + К -Е ьі • кі; і=1

Уп+1 = Уп + К 'ЕЩ • кі;

і=1

I —1

к = /(хп + С  кп;у1 + кп  ЕаН  к1'>1 = 1>->^;

ап+1 = \(уп +1 — уп+1)\

Для определения локальной погрешности менее точного результата и

,п + 1

управления величиной шага интегрирования используется величина а . Вложенный метод Рунге-Кутта р(р) - это схема, в которой метод р -го порядка (“оценщика погрешности”) получается как побочный продукт метода р-го порядка. Порядок аппроксимаций: уп и уп обычно отличается на 1, т.е. р = р +1 или р = р — 1. Оба приближения используют практически одни и те же значения шаговых коэффициентов к1,1 = 1,...,&•, но комбинируют их по-разному. Это позволяет уменьшить количество вычислений функции /(х,у). Для 5- стадийных вложенных методов эта величина составит 5 + с , где с - некоторая константа (обычно с=1 или с = 2).Для сравнения приведем эти оценки для р= 4: накладные расходы для метода Рунге-Кутта с правилом Рунге, а именно число вычислений правой части, равно 11, а для вложенного метода Фельберга того же порядка точности - всего 6.

Метод локальной экстраполяции Ричардсона является обобщением технологии удвоения шага по правилу Рунге [3-5]. Идея этого метода заключается в многократном измельчении шага интегрирования, а также в многократном применении процесса вычисления, названного локальной экстраполяцией. Решение задачи Коши рассматривается при переходе из точки хп в точку хп +1 = хп + Н, Н - базовая длина шага, Н > 0. Выбирается ряд натуральных чисел, такой что: п1 < п2 < ... < пк—1 < пк < ... и, соответственно, последовательность шагов: к1 > к2 > ... > кк—1 > кк > ... , где кI = Н / п^ . Задается опорный числен-

ный метод порядка р и, выполнениям пг- шагов интегрирования длиной кг- вычисляется приближенное решение исходной задачи:

Тг,1 = у^

Выполнив вычисления для ряда последовательных значений I, по рекуррентному соотношению(2) определяют экстраполированные значения Т] ]+1 для произвольных I,] по формуле (4). Этот процесс получил название «локальная полиномиальная экстраполяция»:

Т-  +1 := Т-  | Т,] Т~1]

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

Таблица 1

Экстраполяционная таблица

р р+Ь р+2Ь р +(к — 2 )Ь р + (к — 1)Ь

Т11

Т21 Т22

Т31 Т32 Т33

Тк—1к—1

Тк1 Тк2 Тк3 Тк —1к Ткк

В таблице 1 Т] - приближенное решение задачи Коши, полученное численным методом порядка р + (] — 1)Ь с шагом к-. Величина Тк к соответствует аппроксимации наивысшего порядка, равного 2к, в случае, если вычислены первые к строк экстраполяционной таблицы, а величина Тк—1 к соответствует аппроксимации порядка 2к-2. Для управления шагом интегрирования естественно использовать выражение

||Тк—1,к — Тк,к\\.

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

2.Анализ вычислительной сложности последовательных методов решения СОДУ с постоянными коэффициентами. Исследование эффективности предложенных альтернативных способов оценки локальной погрешности проводилось на следующем множестве методов одного и того же порядка р:

1 - явный метод Рунге-Кутта и правило Рунге;

2 - экспоненциальный метод и правило Рунге;

3 - вложенные методы Рунге-Кутта;

4 - вложенные методы на основе экспоненциальных;

5 - метод Грэгга-Булирша-Штера;

6 - экспоненциальный метод с локальной экстраполяцией.

Экспоненциальный метод относится к специальным методам численного интегрирования линейных СОДУ с постоянными коэффициентами, основанным на точном представлении решения в аналитической форме и вычислении матричной экспоненты. Точным решением задачи Коши вида

Г у' = А • у;

1

[у(хс) = уо

является матричная экспонента

у( хо + к) = еАку о, еАк =ЕЕ -(кА) к

'0 ' ^ о > ^71

к=0 к!

Приближенное решение (6) можно построить, аппроксимировав матричную экспоненту отрезком ряда Тейлора:

еАк ар

к=0 к!

причем (7) определяет численный метод решения (6) порядка р [6].

Наиболее эффективным последовательным методом, реализующим технологию локальной полиномиальной экстраполяции, считается алгоритм Грэгга-Булирша-Штера (ГБШ), базирующийся на модифицированном методе средней точки [3]:

у1 = уо + кк • /уо);

у,+1 = у—1+Щ  /Сх, у У=\...пк —1;

у(х+щ = £ [ уп1 + упк—1 + к • /(х+H, ущ) К = Н / п.

Наиболее трудоемким из рассматриваемых алгоритмов является метод ГБШ, число арифметических операций для него имеет порядок 0([р2 + р + 2]т2), где т - размерность системы уравнений, ар- порядок метода. Наименее трудоемким является метод вложенных форм Рунге-Кутта: 0([2р + 2] • т ) и вложенный метод на основе экспоненты: 0([2р + 6]т ).

3. Анализ эффективности параллельных алгоритмов решения СОДУ с постоянными коэффициентами с учетом локальной погрешности. Рассмотрим отображение алгоритмических схем приведенных методов на многопроцессорные вычислительные системы 81МБ-структуры с распределенной памятью. Конфигурацию системы считаем фиксированной: число процессорных элементов и схема их соединения не изменяются в процессе счета. Каждый процессор может выполнить любую арифметическую операцию за один такт, временные затраты, связанные с обращением к памяти, отсутствуют.

Параллельная реализация перечисленных методов требует распараллеливания следующих базовых операций: матричное и векторное умножение и сложение, а также умножение вектора и матрицы на скаляр. Наиболее оптимальное топологическое решение - это квадратная сетка т х т или ее замкнутый эквивалент - тор. На такой топологической схеме достаточно эффективно выполняются матричные операции. Для простоты изложения рассмотрим случай, когда количество процессорных элементов в строке или столбце матрицы совпадает с размерностью задачи т. Вычисление матричного умножения и умножения матрицы на вектор может быть выполнено по систолическому алгоритму, который является наиболее эффективным для 8ГМЭ-систем [7]. Процесс выполнения систолического умножения матриц состоит из предварительного косого сдвига левого сомножителя по строкам, влево и правого - вверх, по столбцам. Затем пошагово выполняются: т умножений элементов, (т-1) одиночный сдвиг и (т-1) сложение:

тЦа = т!ум +

где !ум , !сл , С - время выполнения одиночных операций умножения, сложения и

сдвига. Вычисление систолического умножения матрицы на вектор на базе алгоритма сдваивания подробно описано в [7] и требует следующего времени выполнения:

ТаШ = !ум +

где т = 21 — 1 и I = \log2m\.

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

Наименее трудоемкими оказались экспоненциальный метод с оценкой погрешности по правилу Рунге, для которого число арифметических операций имеет порядок: 0(3  Log2m) , и вложенный-экспоненциальный метод: 0(2  Log2m) . Худшие показатели потенциального параллелизма у метода Рунге-Кутта с правилом Рунге и метода ГБШ.

Общепризнанным приемом повышения эффективности параллельных вычислений служит сокращение межпроцессорных обменов данными. Уменьшение времени выполнения обменных операций возможно как за счет повышения пропускной способности линков, так и путем целенаправленного программирования решаемых задач, ориентированного на повышение отношения времени обработки к времени передачи данных в ВС. Эта задача тем более актуальна в связи с тем, что скорость обработки данных в микропроцессорах растет существенно быстрее увеличения пропускной способности интерфейсов микропроцессоров [2]. Оценки трудоемкости обменных операций для исследуемых алгоритмов приведены в таблице 2. Наибольшей коммуникационной сложностью обладает метод ГБШ, наименьшей - вложенный экспоненциальный и экспоненциальный на базе локальной экстраполяции.

Таблица 2

Характеристика алгоритмов по трудоемкости обменных операций

Параллельный алгоритм Количество обменных операций

1.Метод Рунге-Кутта и правило Рунге

2. Экспоненциальный метод и правило Рунге 9

3.Вложенные методы Рунге-Кутта 3

4.Вложенный экспоненциальный метод 6

5.Метод ГБШ 3(т -1)(р2 + 2р + 4)/4

6. Экспоненциальный метод с локальной экстраполяцией Ричардсона

Анализ зависимости отношения времени обменов ко времени выполнения арифметических операций приведен на рис. 1.

В целом, характеристики реального параллелизма:

1) лучшие методы:

- 6-экспонента и локальная экстраполяция : 0(0.1т) ;

- 4- вложенная экспонента: 0(0.6т);

2) худший метод - метод ГБШ: 0(0.3  [р2 + 2р + 4]т).

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

Кр(т) = Ті(т)/Тр(т) ,

т.е. к отношению времени решения на скалярной ЭВМ: Т\(т) - к времени выполнения параллельного алгоритма:Тр(т) . Величина т - параметр, количество входных данных задачи, в данном конкретном случае порядок СОДУ.

♦ 1

2

* 3

4

♦ 5

- и— 6

Рис. 1. График зависимости отношения обменных операций к арифметическим

от размерности задачи

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

Кр I

Рис. 2. Зависимость коэффициентов ускорения от порядка системы с учетом

операций обмена

Определение характеристик параллелизма осуществлялось с помощью пакета Mathematica (Wolfram Research Inc.), численный эксперимент проводился на базе тестов для СОДУ, разработанных в НИВЦ МГУ [8] .

Анализ полученных результатов позволяет сделать следующие выводы:

- самым эффективным из рассмотренных последовательных методов является вложенный метод с контролем погрешности на шаге по правилу Рунге;

Z

m

m

- лучшими параллельными методами с точки зрения трудоемкости являются вложенный экспоненциальный метод и экспоненциальный с локальной экстраполяцией;

- преимущества методов на базе локальной экстраполяции проявляются при получении высокоточных решений (10-15 -10-20) и для сложных правых частей.

Заключение. Анализ вычислительной сложности различных технологий определения апостериорной локальной погрешности при решении задачи Коши имел своей целью управление шагом интегрирования в достижении некоторой заданной точности с минимальными вычислительными затратами. Перспективным направлением дальнейших исследований является оценка влияния различных топологий многопроцессорных ВС на характеристики качества параллельных алгоритмов, исследование устойчивости полученных методов.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Rajkumar Buyya. High Performance Cluster Computing. Volume 1: Architectures and Systems. Volume 2: Programming and Applications. Prentice Hall PTR , Pren-tice-Hall Inc., 1999.

2. Воеводин В.В , Воеводин Вл.В. Параллельные вычисления . СПб.: БХВ-Петербург, 2002. 608с.

3. Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений.Нежесткие задачи: Пер. с англ. М.: Мир, 1990. 512с.

4. Bergman S., Rauber T., Runger G. Parallel execution of embedded Runge-Kutta methods.International Journal of Supercomputer Applications, 10(1):62-90, 1996.

5. Houwen P.J., Sommeijer B.P. Parallel ODE solver.// Proceedings of the International Conference on Supercomputing. ACM Press, 1990, p. 71-81.

6. Арушанян О.Б., Залеткин С.Ф. Численное решение обыкновенных дифференциальных уравнений на Фортране. М.: МГУ,1990. 336с.

7. Бройнль Т. Паралельне програмування:Початковий курс: Навч. поабник/ Вступ.слово А.Ройтера;Пер. з шм. В.А.Святного. К.: Вища школа, 1997. 358с.

8. Арушунян О.Б., Залеткин С.Ф., Калиткин Н.Н. Тесты для вычислительного практикума по обыкновенным дифференциальным уравнениям.// Вычислительные методы и программирование, 2002.Т.3. С. 11-19.

А.В. Шматок, З.В. Каляев

СТРУКТУРНО-ФУНКЦИОНАЛЬНЫЙ И СТРУКТУРНОДЕКЛАРАТИВНЫЙ МЕТОДЫ РАЗРАБОТКИ КОМПОНЕНТОВ МАСШТАБИРУЕМЫХ ПАРАЛЛЕЛЬНЫХ ПРОГРАММ

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