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

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

Вестн. Ом. ун-та. 2011. № 4. С. 157-163.

УДК 517.91

В. В. Коробицын, Ю.В. Фролова

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

Предлагается параллельная реализация метода последовательных приближений с использованием технологии МУГОІЛ СИБЛ. Приводятся результаты тестирования алгоритма на видеокарте. Показано влияние значения параметров метода на точность и устойчивость. На примере одной задачи показан способ получения оптимальных значений параметров, обеспечивающих наибольшую скорость вычисления траектории решения системы обыкновенных дифференциальных уравнений.

Ключевые слова: метод Пикара, параллельный алгоритм, технология СИБЛ.

Введение

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

Одной из актуальных задач параллельных вычислений является разработка параллельных алгоритмов численного решения систем обыкновенных дифференциальных уравнений (ОДУ). Однако большинство алгоритмов решения систем ОДУ являются последовательными. Поэтому наиболее распространенный вариант распараллеливания заключается в параллельном вычислении правой части системы [1; 5]. Особенно это выгодно, когда размерность системы велика (более 1000) [3; [4]. Такое распараллеливание принято называть распараллеливанием по пространству (parallelism across space). Однако существует и другой подход, который реализует распараллеливание вдоль траектории решения системы. Траектория вычисляется одновременно в нескольких точках. Такое распараллеливание принято называть распараллеливанием по времени (parallelism across time). Одним из первых алгоритмов, который можно так распараллелить, является алгоритм последовательных приближений (метод Пикара).

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

© В.В. Коробицын, Ю.В. Фролова, 2011

циалированных процессорах НУЮ1Л.

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

1. Метод последовательных приближений

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

В общем виде метод формулируется так. Необходимо расщепить функцию Л^х) из дифференциального уравнения

&х

&

= / (и х) = Ш, х) + /2^, х), х (^) = Х0

так, чтобы легко можно было бы решить уравнение вида Сх/dt=/l(t,Х)+c|(t). Тогда, начиная с первого приближения Хэ(^, на каждой 1-ой итерации решаем уравнение

&х'+1

-— = х1+1(0) + х‘(0), х1+1(0 = х0, (1)

&

получаем последовательность функций х1^), х2(^,..., равномерно сходящуюся к решению х(^.

В простейшем случае выбирают / = 0, / = /, тогда, проинтегрировав (1), получаем формулу

х'+1(t) = хо + I /(5, х' (5))&5.

(2)

Этот вариант метода имеет название

- метод последовательных приближений Пикара. Последовательность х'(0 равномерно сходится к решению х(0 при вы-

полнении условия Липшица для функции //,х) по второму аргументу.

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

При использовании методов Рунге-Кутты невозможно вычислить решение в точке t+h, если неизвестно значение в точке t. Тогда как при вычислениях по методу последовательных приближений при нахождении интеграла значения функции / можно вычислять параллельно для различных точек t, что способствует построению параллельного вычислительного алгоритма.

2. Численная схема метода последовательных приближений

Будем рассматривать задачу Коши для системы из п обыкновенных дифференциальных уравнений первого порядка. Для заданной функции ГК^^К” необходимо найти х:[70; Т] ^К”, удовлетворяющую системе ОДУ

^ 1 (t, х) (3)

&

с заданным начальным условием х(^)=х0.

Предполагаем, что функция 1^, х) удовлетворяет условию Липшица по второму аргументу во всей области определения, тогда решение существует и единственно.

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

следовательность узлов ^, ^, ^,..., tm.1 = Т. Обозначим за х/ приближенное значение функции х(0 в точке tj на итерации '. В качестве начального приближения возьмем х/ = х0, то есть на нулевой итерации функция х°(0 постоянная. На каждом интервале [г); tj+\] решение будем находить с помощью итерационной процедуры (2). В этой формуле присутствует определенный интеграл, который будем вычислять численно.

В зависимости от выбора квадратурной формулы получим методы различной точности. В частности, если выбрать ме-

тод трапеций, то получим следующую итерационную процедуру:

sj=(f(0,xr)+f(tj+i,x-)^ m

xj,i= xj + sj, (5)

x;- = x0, x0! = x0 для всех j = 0,i, ..., m-2, i = i,2,...,K.

Обращаем внимание на то, что в формуле (4) при вычислении sji отсутствуют значения x на i итерации, что дает возможность производить вычисления для i итерации независимо для всех точек tj одновременно, используя только значения на предыдущей итерации i-i. Этот факт дает возможность параллельной реализации вычисления формулы (4) для всех j одновременно. Вычисление формулы (5) реализуется последовательным алгоритмом. Таким образом, вся итерационная процедура (4)-(5) разбивается на параллельную и последовательную части алгоритма. Но поскольку в последовательной части присутствует только одна операция сложения, а в параллельной -трудоемкое вычисление функции f, предложенный алгоритм должен иметь значительное увеличение скорости вычислений, несмотря на то, что потребуется произвести K итераций.

3. Параллельная реализация численного метода на CUDA

Вычислительная производительность современных видеокарт во много раз превосходит производительность центральных процессоров. Архитектура графических процессоров позволяет производить большое количество параллельно исполняемых потоков вычислений на одном GPU. Кроме того, компания NVIDIA выпускает платы Tesla, предназначенные исключительно для решения вычислительных задач, поддерживающих архитектуру CUDA. Существуют кластерные решения, объединяющие серверы на основе х8б архитектуры и вычислительные узлы c Tesla, так называемые гибридные системы. Для реализации метода последовательных приближений была выбрана технология NVIDIA CUDA.

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

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

Конфигурация сетки ядра задается параметрами: r - количество вычисляемых траекторий решения, n - размерность системы ОДУ, m - количество шагов по времени в схеме метода. Входные параметры функции: массив х, содержащий векторы начальных данных (r x n элементов), t0 - начальный момент времени, h - величина шага по времени (|h|>0), K - количество итераций приближения (K > 2), stp - количество больших шагов длиной H = (m - 1)h. Выходной параметр -массив у, содержащий векторы решения в момент времени t0 + stpH.

Для обеспечения высокой скорости обмена между потоками эффективно использовать разделяемую память GPU (shared memory). Эта память имеет высокую скорость доступа, но ограничена в объеме. В предлагаемой реализации используется 2 массива (р и s), размещенные в разделяемой памяти. Массив р используется для хранения значений xj, массив s

- для хранения значений правых частей f(t, x;) (j = 0,...,m - 1). Поскольку каждый элемент Xj имеет размерность задачи n, то размерность массива получается m x n. Для хранения чисел выбран типа float.

Функция ядра odesa реализует следующие этапы:

1) определение индексов потока и занесение начальных значений х0 во все ячейки массива p;

2) вычисление stp шагов длины H по схеме метода;

3) копирование значений найденного решения в массив y.

На каждом вычисляемом шаге (этап 2) производится K итераций, осуществляющих действия:

• вычисление функции правой части системы fk(tj, х/"1) и сохранение ее значений в массиве s;

• вычисление Ху+і по формулам (4)-(5) и сохранение этих значений в массиве р;

• копирование вычисленных значений хКт-1 во все ячейки массива р после завершения К итераций.

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

Вызов функции ядра ойеъа осуществляется из функции ойе$а_$ої\е. Данная функция является интерфейсом для между СРи и ОРи и осуществляет следующие действия:

• выделяет память на устройстве;

• задает конфигурацию ядра;

• копирует начальные данные задачи на устройство;

• запускает ядро;

• копирует результат вычислений с устройства в оперативную память;

• освобождает память на устройстве.

Функция odesa_solve сохраняет результат вычислений в массиве х.

4. Теоретическая оценка

показателей эффективности параллельной реализации

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

Для вычисления г траекторий системы ОДУ размерности п для т шагов интегрирования методом Эйлера потребуется время

Т1(г,п,т) = г п т (/+ ґа + іД где ґ/ - среднее время для вычисления одной компоненты вектор-функции правой части системы, іа и - время выполнения операций сложения и умножения.

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

Т 1(г,п,т) = К г п т (/+ 2 ґа + іД

При параллельной реализации на р процессорах имеем

Т р(г,п,т) = К г п т ґ/+2 ґа + ґр + іс) / р, где їо - среднее время синхронизации потоков. Получаем теоретическую оценку ускорения параллельного алгоритма

Sp (r, n, m) =

Tj*( r, n, m) Tp( r, n, m)

= p

i —

t f + 2t +1 +1

f а ц а

Если ^ мало по сравнению с / то ускорение будет стремиться к р. Если сравнить производительность параллельного алгоритма с последовательным методом Эйлера, то получим ускорение

З ( ) т (г,”,т) p(t/+а+

Зр (г, ”, т) = ^Л---------------------*---

Tp (r,n, m)

K (tf + 2ta+ tu+ ta)

/ ■ а ' -р

В этом случае, Бр<р/К, а эффективность применения метода последовательных приближений при параллельной реализации на р процессорах будет

I p (r, n, m) =

Sp (r, n, m) i

K

<

P

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

При реализации алгоритма на различных аппаратных средствах необходимо учитывать не только количество процессоров, но и их производительность (количество выполняемых арифметических операций с плавающей точкой в секунду, flops). Обозначим производительность последовательного процессора v1, а производительность параллельного процессора vp. Тогда теоретическое ускорение параллельной реализации метода последовательных приближений на параллельном процессоре относительно метода Эйлера в последовательном исполнении определяется формулой

v

Sp=kV- . (6)

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

5. Численные эксперименты

Исследование численного алгоритма проводилось на нелинейной модели динамики численности двух взаимодействующих популяций (упрощенная модель Лот-

ки-Вольтерра). Модель описывается системой из двух обыкновенных дифференциальных уравнений:

dx- 2

—1 = ax, - ex, x2 - ax, , dt 1 1 2 1

dx

—2 = уx2 + ^x, x2 - bx\. dt

(7)

x 2

x -

Рис. 1. Траектория решения системы (7) на фазовой плоскости

В модели использовались следующие параметры: a = 0,25, в = 0,01, у = -0,3,

a = 0,001, b = 0,002. Численное решение вычислялось для траектории с начальными данными: x1(0) = 50, x2(0) = 30 на интервале времени te[0; 100]. Таким образом, были заданы параметры задачи: n = 2, t0 = 0, T = 100, х[0] = 50, х[1] = 30. Траектория на фазовой плоскости (рис. 1) показывает поведение решения задачи. Для оценки производительности метода вычислялось одновременно 8192 траектории. Исследование алгоритма проводилось для следующих значений параметров метода: те{4; 8; 16; 32; 64; 128; 256}, Ke{2; 3; 4; 5; 6; 7; 8; 9}, he[0.00001; 1.0].

Для вычислений мы использовали тип данных float поскольку точность выбранной схемы интегрирования (метод трапеций) невысока. Поэтому использование типа double потребовало бы больше вычислительных затрат с достижением той же точности.

Сначала мы оценили погрешность численного решения при разных заданных параметрах метода. Графики зависимостей (рис. 2) глобальной погрешности E найденного численного решения от выбранного шага интегрирования h при различных значениях параметра m (от 4 до 256) и фиксированном K = 3 показывают поведение погрешности метода. Глобаль-

ная погрешность E оценивалась как максимум отклонений численного решения Xj от точного x(tj) для всех j = 0,1,.,N: E = max_j |х;- - x(tj)|, здесь N - общее количество вычисленных точек траектории решения. В эксперименте значения точного решения мы вычисляли с помощью метода Рунге-Кутты-Мерсона 5 порядка с шагом, достаточным для получения решения с точностью, требуемой для типа данных float. Далее в тексте под погрешностью мы будем подразумевать глобальную погрешность численного решения, вычисленную именно таким способом.

Рис. 2. Графики зависимости погрешности численного решения Е от выбранного шага интегрирования Л для метода Эйлера (1) и метода последовательных приближений при заданных значениях параметров метода К = 3 и различных т: т = 4 (2), т = 16 (3), т = 64 (4), т = 256 (5)

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

Приведенные графики демонстрируют поведение погрешности численного решения вблизи границы области сходимости. Когда шаг интегрирования превышает граничное значение области сходимости, происходит быстрый рост погрешности. Обозначим эту величину за h . Так, для т =16 значение на границе области сходимости h «0,12 (-^ h ~ 0,9), для т = 64 получили h « 0,016 (-^ h « 1,8), для т = 256 - к* « 0,002 (-^ к* « 2,7). Когда h < к* погрешность зависит линейно от h. На исследуемой задаче мы получили зависимость Е = 0,05k в области сходимости и

Е = 0,02т2к2 - 0001тк + 0,001 вне области сходимости. Из графиков видно, что в области сходимости метода погрешность решения немного меньше, чем у метода Эйлера.

Рис. 3. Графики зависимости погрешности Е численного решения от выбранного шага интегрирования Л для метода последовательных приближений при заданных значениях параметров метода т = 256 и различных К:

К = 2 (1), К = 3 (2), К = 4 (3), К = 6 (4), К = 9 (5)

Далее мы оценили влияние параметра К на область сходимости метода. Графики (рис. 3) зависимости погрешности Е от длины шага интегрирования при фиксированном параметре т = 256 и различных значениях параметра К (от 2 до 9) показывают влияние параметра К на сходимость метода. Увеличение параметра К приводит к расширению области сходимости метода. Так, для К = 3 граничное значение к = 0,002, для К = 4 получаем к = 0,004, для К = 6 имеем к = 0,016, для К = 9 достигаем к = 0,025. Дальнейшее увеличение параметра К дает незначительное увеличение области сходимости метода по сравнению с ростом объема вычислений.

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

tc(Е") = тттк .

к ,Ек < Е*

Аппаратное обеспечение для тестирования включало компьютер с центральным процессором Intel Core2Quad Q6600 2,40 ГГц с пиковой производительностью 38,4 Gflops (по 9,6 Gflops на каждое ядро) и видеокартой с графическим процессором NVIDIA GeForce 8800GT с пиковой производительностью 351,7 GFlops. В данной конфигурации оборудования получаем теоретическую оценку ускорения Sp = 18,3, вычисленную согласно формуле (6) при K = 2, vp = 351,7, v1 = 9,6.

lg tc

Рис. 4. Графики зависимости наименьшего времени вычисления (о 8192 траекторий от выбранного параметра метода т при разных заданных предельных погрешностях Е численного решения:

Е = 0,01 (1), Е = 0,001 (2), Е = 0,0001 (3)

Графики (рис. 4) зависимости минимального времени (с от параметра т при

Г*

различных порогах точности Е показывают нелинейную зависимость 4 от т. Время измерялось для вычисления 8192 траекторий задачи (7). Оптимальные значения, полученные в результате проведенного эксперимента, приведены в таблице. Здесь же для сравнения приведено время Т1 решения той же задачи с помощью метода Эйлера на центральном процессоре. Реальный показатель ускорения = Т1Пс указывает максимально достигнутое ускорение для заданной точности. Полученные результаты показывают 15-тикратное ускорение при теоретической оценке в 18,3, что является весьма приемлемым результатом.

Таблица 1 Оптимальные значения параметров метода

E m K h E tc Ti S

0,01 8 2 0,5 0,008 0,016 0,265 15,9

0,001 8 2 0,06 0,0007 0,140 2,20 15,7

0,0001 8 2 0,008 0,00008 1,203 17,6 14,6

Выводы и заключение

Проведенный эксперимент с параллельной реализацией метода последовательных приближений показал: 1) метод пригоден для распараллеливания на графическом процессоре; 2) точность приведенной реализации метода невысока (глобальная погрешность линейно зависит от шага интегрирования к); 3) область сходимости метода уменьшается с увеличением параметра т; 4) увеличение параметра К расширяет область сходимости метода; 5) оптимальные значения параметров метода т = 8, К = 2, шаг к опреде-

1-*

ляется исходя из заданной точности Е ; 6) достигнутое ускорение решения задачи относительно метода Эйлера в последовательном исполнении составило 15 раз.

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

ЛИТЕРАТУРА

[1] Новиков Е. А., Каменщиков Л. П. Реализация полуявного (4,2)-метода решения жестких систем на параллельной ЭВМ // Вычислительные технологии. 2004. Т. 9. Вестник КазНУ им. аль-Фараби. Сер. Математика, механика, информатика. № 3 (42). (Совм. выпуск. Ч. 1). С. 235241.

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

[3] Abbas S.H. A large parallel block predictor-corrector method for initial value problems // Chaos, Solitons and Fractals. 2008. Vol. 38. № 1. P. 263-269.

[4] Ackermann J., Baecher P., Franzel T., Goese-le M., Hamacher K. Massively-Parallel Simulation of Biochemical Systems // Proceedings of Massively Parallel Computational Biology on GPUs, Lubeck, Germany, September 29, 2009.

[5] Amodio P., Brugnano L. Parallel solution in time of ODEs: some achievements and perspectives // Applied Numerical Mathematics. 2009. Vol. 59, № 3-4. P. 424-435.

[6] Gear C.W., Xuhai X. Parallelism across time in ODEs // Applied Numerical Mathematics. 1993. Vol. 11. № 1-3. P. 45-68.

[7] Коробицын В. В. Программная реализация

метода последовательных приближений на CUDA. URL: users.univer.omsk.su/~korobits/

cuda_app.