Алгоритм реализации параллельных вычислений по формулам численного интегрирования Рунге-Кутта Текст научной статьи по специальности «Компьютерные и информационные науки»

УДК 658.512:004.42,658.512:519.87

К.Г. Жуков

алгоритм реализации параллельных вычислении по формулам численного интегрирования рунге-кутта

Вычислительные структуры для потоковых задач

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

В общем случае потоковую задачу можно сформулировать следующим образом. Предположим, что существует упорядоченное множество векторных данных < , , ..., >, (1 = 1, 2, ..., N), каждый элемент которого должен быть обработан по фиксированному алгоритму, представляемому в виде графа G(Q, X) (рис. 1).

Каждой вершине qj е Q графа G(Q, X) поставлена в соответствие операция О, входящая во множество допустимых операций O. Любая дуга (qj, qj+1) е X определяет, что результат операции

0. является входным операндом для операции

01, поставленной в соответствие вершине q1.

Задача потоковой обработки состоит в преобразовании потока векторов входных данных Д (1 = 1, 2, ..., N) в поток векторов выходных данных B¡ <Ь,Ь2, ..., Ь'ь >,(1 = 1, 2, ..., Щ в соответствии с графом - алгоритма G(Q, X). Решение поставленной задачи можно выполнить на после-

Рис. 2. ЭВМ с фон-неймановской архитектурой

довательной ЭВМ с фон-неймановской архитектурой (рис. 2).

Процессор Р необходимо запрограммировать на последовательное выполнение операций 0(1 = 1, 2, ..., М), соответствующих вершинам графа G(Q, X), и затем обеспечить ввод последовательности векторов (¡' = 1, 2, ..., N).

Время обработки N векторов входных данных определяется выражением М

Греш= N I ¡(0 )Т, (1)

1 = 1

где M=\Q\ - число вершин на графе-алгоритма G(Q, X); ¡(О¡) - число тактов работы ЭВМ при выполнении операции О., соответствующей вершине графа-алгоритма G(Q, X); т - продолжительность такта.

Критерием оценки выполнения задачи является условие Треш<Трв. Здесь Трв - заданное время на обработку всех векторов входных данных. Если условие выполняется, то такой вариант реализации следует признать эффективным с точки зрения простоты реализации. Однако при достаточно малых значениях Т__ условие Т < Т__ не

РВ ^ реш РВ

будет выполняться.

Рис. 1. Граф алгоритма

Рис. 3. Параллельный способ обработки векторов данных

Очевидным путем сокращения времени обработки Ореш входных массивов является распараллеливание процесса обработки. В простейшем случае такой способ подразумевает наличие С процессоров Р. (г = 1, 2, ..., С), каждый из которых может работать независимо от других процессоров. Разумеется, что если каждый из процессоров Р. (г = 1, 2, ..., С) запрограммирован на реализацию графа-алгоритма О(0, X), а множество входных векторов Д (г = 1, 2, ..., N) разбито на ШС непересекающихся подмножеств, то каждое из подмножеств может быть обработано на процессоре Р. независимо, т. е. параллельно с другими подмножествами (рис. 3).

В этом случае время решения будет составлять:

М

Треш = N / С I /(О. )т. (2)

г=1

По сравнению с (1) Треш сокращается в С раз.

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

рывного поступления вектора входных данных (г = 1, 2, ..., N) на первую ступень конвейера, состоящего из а ступеней (процессоров). Основное преимущество конвейерного способа по сравнению с параллельным - сокращение числа входных и выходных каналов - для многих задач является несущественным. Для таких задач число операций М на графе значительно больше числа входных векторов N.

Естественный путь увеличения скорости обработки информации в конвейере - это распараллеливание обработки на каждой его ступени (рис. 4).

Время обработки потока из N векторов на этом конвейере будет составлять

T = (N + H -1)* AT,

реш 4 / '

(3)

причем для данного случая

AT = mах(т;еШ) = max(max f (Oj)* т) =

i=1...H F i=1...H j=1...M J

= max f (O )* т,

i=1...M

где max f (O) - наиболее длительная опера-

i=1...M

ция, входящая в состав вершин графа G(Q, X).

Рис. 4. Представление г-й ступени конвейера в виде М. параллельно работающих процессоров

При большом значении N, а также полагая, что М

М1 «—, выражение (3) можно представить в Н

приближенном виде:

I/ (О):

!=1

Т = N■

реш Н * М / Н

N М

=Н I / (0j у

Н j=1

т.

(4)

Таким образом, время решения по сравнению с обычным конвейером сокращается приблизительно в — раз, где М - число вершин в графе Н

G(Q, X); Н - число ступеней конвейера.

Такая схема будет эффективна только при условии, что все вершины подграфа G(Q¡, X) (1 = 1, 2, ..., Н) информационно независимы и могут быть реализованы параллельно (одновременно).

В общем случае операционные вершины подграфов G(Q¡, X) являются информационно зависимыми, т. е. некоторая вершина qj е Q¡ не может быть реализована до тех пор, пока не получен результат реализации вершины qj_1 е Qi. Поэтому максимальный темп конвейерной обработки и, как следствие, Треш обработки всего множества входных данных можно обеспечить, если структура связей между ступенями конвейера будет полностью адекватна топологии связей вершин в графе-алгоритма G(Q,X).

Такой способ организации конвейерных вычислений называется структурным. Максимального темпа обработки можно достичь, если каждой операционной вершине q¡ (¡' = 1, 2, ..., М) графа-алгоритма поставить в соответствие свой процессорный элемент Р. ^ = 1, 2, ..., М) и обеспечить связи между Р. согласно топологии информационных дуг графа G(Q, X) (рис. 5).

Если на вход такой вычислительной структуры последовательно подавать вектора

Б. (¡' = 1, 2, ..., N входного множества, то время их обработки составит Треш = ^ + Нтах _ 1) * ДТ, где Нтах - число вершин на критическом пути в графе G(Q, X). Критическим путем являет-

М

ся тот путь, для которого значение I /(0j) * т

j=l

максимально; О^ = 1, 2, ..., К) - операции, принадлежащие вершинам критического пути; ДТ = тах( / (О )) * т - время выполнения наибо-

j=1...Н^ ■>"

лее длительной операции, принадлежащей вершинам критического пути. При большом значении N и учете того, что

I / О)

ДТ = тах / (О ) « N

¡=1,2,..., Нт 1

j=1

М

получаем

N М

Т =— I / (О,. )*х. реш М^^

(5)

!=1

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

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

Рис. 5. Структурные конвейерные вычисления

Основные задачи статьи

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

математическому описанию разработанного алгоритма параллельных вычислений по формулам Рунге-Кутта четвертого порядка точности;

созданию в среде LрbVIEW виртуального параллельного решателя;

сравнению структуры вычислителя (решателя) со структурами решателей, используемых в инструментальной среде моделирования LрbVIEW;

тестированию параллельного решателя;

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

Математическое описание разработанного алгоритма

В настоящее время решение обыкновенных дифференциальных уравнений (ОДУ) на ЭВМ с фон-неймановской структурой выполняется методом Рунге-Кутта четвертого порядка. Метод численного интегрирования (ОДУ) разработан немецким математиком Рунге и модифицирован Кутта в начале ХХ в. [2].

Рекуррентные формулы для уравнения

йу / йг с /(У, г), у(0) с у (6)

имеют следующий вид:

= /(у,, г,); (7)

к2 = /(у, + Ъ/2* к„ г, + Ъ/2); (8)

¿3 = /(У + Ъ/2* кг, г, + Ъ/2); (9)

К = /(У, + ъ * ¿3,г1 + Ъ); (10)

У,+1 = У, + Ъ / 6*(£, + 2* ¿2 + 2* ¿3 + £4) с с у, + Ъ /6% + Ъ/3*£2 + Ъ/3*£3 + Ъ/6*£4:

На каждом ,-м шаге интегрирования необходимо производить последовательные вычисления по формулам (7)—(11). Последовательная структура алгоритма вычислений исключает возможность организации параллельного (одновременного) вычислительного процесса. Выполним преобразование выражения (11) путем подстановки ¿1 в правую часть (8); вычисленного значения ¿2 в правую часть (9) и т. д.

Для уравнения (6) с /(у, г) = а * у + Ь * х(г) преобразованное выражение (11) будет иметь следующий вид:

(11)

(13)

у,+1 с у1 + Ъ / 6*(а * у1 + Ь * х1) + +2* Ъ /3*[а * у1 + Ь *( х + Ъ/2)] + +Ъ / 6 *[а * у, + Ь * (х, + Ъ)] + +а*Ъ2 /6* (а* у1 + Ь * х) + (12) +2* а * Ъ2 / 6*[а * у, + Ь *(х, + Ъ/2)] +

+а * Ъ3 /12 * (а * у1 + Ь * х,) + +а2 *Ъъ /12*[а* у1 + Ь *(х, + Ъ/2)] + +а3 *Ъ4 /24 * (а* у1 +Ь * х,.).

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

у,+1 с у, + (а1* у, + Ь1* х,) +

+[а2 * у, + Ь2*(х, + Ъ /2)] + +Ц * у, + Ьх* (х, + Ъ)] + (а3 * у, + Ь3 * х,) + +[а4* у, + Ь4 * (х, + Ъ /2)] +

+(а5 * у1 + Ь5 * х,) + [а5 * у, + Ь5 * (х, + Ъ /2)] +

+(а6* у, + Ь6* х,),

где коэффициентам а1,Ь1, а2, Ь2, ..., а6, Ь6 присвоены значения:

а1 с а * Ъ /6; Ь1 с Ь * Ъ /6;

а2 с а*2*Ъ/3; Ь2 сЬ *2*Ъ/3; а3 с а2* Ъ2/6; Ь3 с Ь * а * Ъ2/6; а4 с 2*а2 *Ъ2 /6; Ь4 сЬ *2*а*Ъ2; а5 с а3* Ъ3 /12; Ь5 с Ь * а2* Ъ3/12; а6 с а4* Ъ4 /12; Ь6 с Ь * а3* Ъ4/12.

Выполним сравнение двух вариантов реализации формул Рунге-Кутта по числу выполняемых операций умножения и сложения на ,-м шаге интегрирования. Для классического способа (см. (7)—(11)) необходимо 15 операций умножения и 11 операций сложения. Новый алгоритм (см. (13)) требует выполнения 16 операций умножения и сложения. Можно заметить незначительный рост в числе операций (6 операций).

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

Создание виртуального решателя

Реализацию решателя ОДУ выполним в ин-

Рис. 6. Блок-диаграмма параллельного векторного решателя

струментальной среде разработки виртуальных приборов (ВП) LabVIEW 2010 фирмы National Instruments. Выбор среды LabVIEW обусловлен несколькими важными ее свойствами:

графическим способом программирования; потоковым принципом вычислений; возможностью выполнения векторных операций;

наличием встроенного решателя ОДУ методом Рунге—Кутта четвертого порядка.

Тестирование решателя проведем на примере ОДУ

dy/dt = a*y + b*t c y(0) = 1 и a = b = 1. Блок-диаграмма виртуального решателя представлена на рис. 6. В левой части блок-диаграммы производится формирование шести векторов коэффициентов [a1 bj, [a2 b2L [a6 b6] и трех вект°р°в |y0 x0L [y0 x0+h/2], [y0 x0+h] с компонентами, определяющими начальное значение y0 зависимой переменной y и значения возмущающей функции (x(t) = t)x0, х0 + h/2, х0 + h в начале, середине и в конце промежутка интегрирования.

В теле цикла размещены пять процессоров. Четыре процессора (Процессор 1, ..., Процессор 4) построены на базе восьми векторных умножителей, восьми векторных сумматоров и трех скалярных сумматоров. На выходах процессоров формируются значения приращений решения y которые суммируются со значением y на скалярном сумматоре, вычисляющем значение y.+p (см. (13)). Последний сумматор входит в со-

став Процессора 5, дополнительно формирующего новые векторы [y1 x1],[y1 x1 + h/2],[y1 x1 + h].

Основной цикл реализован на основе структуры LabVIEW For Loop. Для выполнения итерационных вычислений в структуре предусмотрены так называемые сдвиговые регистры (Shift Registers). Сдвиговые регистры обеспечивают хранение и перезапись значений переменных, например, при реализации присваивания 1 после окончания i-й итерации. Регистры допустимо применять и при работе с векторами. В блок-диаграмме используется семь сдвиговых регистров.

Тестирование векторного решателя

Для оценки правильности выражений (13) и проверки работоспособности виртуального векторного решателя в состав блок-диаграммы включен функциональный блок ODE Runge—Kut-ta 4-th Order.vi, входящий в состав библиотеки математических функций NI_Gmath.lvlib. Блок-диаграмма стандартного решателя представлена на рис. 7.

На блок-диаграмме отчетливо отражена реализация последовательности выражений (7)—(11). С помощью четырех функциональных блоков R-K f (x,y) вычисляются значения кх, к2,k3 ,кА. Узел Formula Node реализует часть выражения (11) (без выполнения умножения на значение шага интегрирования h). Вычисления завершают блок умножения (h*phi) и сумматор (л = х0 + +h*phi). Данные представлены в основном в век-

Рис. 7. Блок-диаграмма стандартного решателя LabVIEW

торном виде. Векторная реализация обеспечивает решение систем дифференциальных уравнений. Значения к1, к2, к3, к4 в каждом из четырех блоков R-K /(х, у) вычисляются для всех уравнений системы «параллельно». Результаты интегрирования тестового уравнения представлены на лицевой панели виртуального векторного решателя (рис. 8).

В окне Waveform Graph изображены совпадающие графики изменения погрешности двух вариантов построения решателей (Error Y Error Y2). Точные значения вычисляются по аналитическому решению тестового уравнения y(t) = 2* e' -1 -1. Максимальные значения погрешностей для шага интегрирования h = 0,01 на всем отрезке интегрирования

Рис. 8. Лицевая панель векторного решателя

0 < = t < = 5 не превосходят величины 1,226516701535729E-7.

Элементная база для реализации параллельного решателя

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

Для аппаратной реализации решателя наиболее подходят кристалы ПЛИС (FPGA) Virtex 5 фирмы Xilinx и Cyclone 4 фирмы Altera. ПЛИС предназначены для создания устройств цифровой обработки данных (DSP) с поддержкой выполнения векторных операций в форме с плавающей точкой (Floating Point IP Cores).

СПИСОКЛ

1. Каляев, И.А. Реконфигурируемые мульти-конвейерные вычислительные структуры [Текст] / И.А. Каляев, И.И. Левин, Е.А. Семерников, [и др.]; 2-е изд. —Ростов-на-Дону: ЮНЦ РАН, 2009. —С. 20—32.

Доказана правильность алгоритма параллельных вычислений по формулам численного интегрирования ОДУ методом Рунге—Кутта.

Погрешности решений тестового уравнения на предложенном решателе и встроенном решателе LabVIEW совпадают.

Алгоритм применим к другим и-этапным методам Рунге—Кутта

Для эффективной аппаратной реализации предложенного решателя необходимо применять ПЛИС (FPGA) Virtex 5 или Cyclone 4 с поддержкой Floating Point IP Cores.

Время решения ОДУ на предложенном решателе в четыре раза меньше, чем на встроенном решателе LabVIEW (и Matlab).

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

ГЕРАТУРЫ

2. Эдвардс, Ч.Г. Дифференциальные уравнения и краевые задачи: моделирование и вычисление помощью Mathematica, Maple и Matlab [Текст] / Ч.Г. Эдвардс, Д.Э. Пенни; Пер. с англ.; 3-е изд. —М.: ИД «Вильямс», 2009. —С. 203—217.

УДК 53.096,536.58

В.Н. Козлов С.В. Хлопин

разностные схемы для си Нтеза упр а в лениЯ нелинейными теплопроводяЩИМИ объектАМИ

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

( du

И U 1 = И 2

dx ax3 [ и" 4(u) (дЯ ax5(u)

+Ay 2

(

dy

Ay 3

д

Ay4(u) | dy Ay5(u)

(1)

+А

dZ az31 Az4(u) [dZ Az5(u)

+ f (x, y, z, t),

и(0, x, y, г) = и0(^ y, г), и ^у,г^ = и3 .

Это уравнение позволяет учесть изменения теплофизических свойств среды в виде кусочно-квадратичных функций (операторов) или их аппроксимаций (кусочно-линейных функций). Эти функции (операторы) позволяют учесть зависимость свойств теплопроводности от совокупности переменных: от первых производных температуры по времени у = фД)); вторых производных температуры от сложной координаты у = ф2 (г (?)); первых производных температуры от совокупности координат оператором у = ф3 ());