Южно-Уральский государственный университет В работе описывается параллельный алгоритм решения нестационарных задач линейного программирования большой размерности, ориентированный на кластерные вычислительные системы. В основе алгоритма, получившего название «следящий», лежат фейеровские отображения. Алгоритм отслеживает изменения исходных данных и вносит корректировки в вычислительный процесс. При этом задача разбивается на большое количество подзадач, которые могут решаться независимо без обменов данными. Приводятся диаграммы деятельности UML, описывающие параллельный следящий алгоритм.
Определим отображение \varphi : \BbbR n
\BbbR n следующим образом: (1) (2) Пусть M – многогранник, задаваемый ограничениями задачи линейного программирования (1). Такой многогранник всегда является выпуклым. Известно [8], что \varphi будет однозначным непрерывным M -фейеровским отображением для любой системы положительных коэффициентов \{ \alpha j > 0\} (j = 1, . . . , m), таких, что \sum \alpha j = 1, и коэффициентов релаксации 0 < \lambda j < 2. Полагая в формуле (2) \lambda j = \lambda
и \alpha j = 1/m (j = 1, . . . , m), получаем формулу
(3) \| aj\|
2 s r\ightaow\infty которая будет использоваться в следящем алгоритме.
Обозначим \varphis(x) = \varphi . . . \varphi(x).
u\nderbac{} s
u\nderbac{} +\infty \{ \varphis(x0)\} s=0 \rightaow x=\ \in M. Под фейеровским процессом , порождаемым отображением \varphi при произвольном начальном приближении x0 \in \BbbR n, будем понимать последовательность \{ \varphis(x0)\} s=0 занный фейеровский процесс сходится к точке, принадлежащей множеству M: +\infty . Известно, что укаБудем кратко обозначать это следующим образом: lim \varphis(x0) = x\=.
Под \varphi-проектированием (псевдопроектированием ) точки x \in \BbbR n понимается отображение \piM \varphi : \BbbR n r\ightaow
M [9], задаваемое соотношением \pi M\varphi (x) = lim \varphis(x).
на многогранник M s r\ightaow\infty 3. Следящий алгоритм
В общем виде следящий алгоритм может быть описан следующим образом. Все пространство \BbbR n делится гиперкубической сеткой на ячейки с переменной длиной грани. Размер ячейки может динамически меняться в зависимости от скорости изменения исходных данных: чем быстрее меняются данные, тем больше становится ячейка, чтобы успевать следить. Алгоритм постоянно находит (отслеживает) ячейку наименьшего размера, на которой достигается максимум целевой функции. Будем такую ячейку называть целевой. В каждой итерации алгоритма просчитываются ячейки, входящие в некоторую следящую область вокруг целевой ячейки. В простейшем случае следящая область может быть кубической формы со следящей ячейкой в центре. В общем случае следящая область может быть произвольной выпуклой фигурой. При этом ячейки всегда имеют кубическую форму. Целевая точка – это некоторая точка, находящаяся вне следящей области. Ее координаты могут динамически вычисляться по коэффициентам целевой функции. Например, при целевой функции Q (x) = x1 + 2x2 в качестве целевой точки можно взять точку (T, 2T ), где T – достаточно большое положительное число. Нулевой вершиной кубической области будем называть вершину куба, находящуюся ближе всего к началу координат. Нулевой вершиной кубической ячейки будем называть вершину ячейки, находящуюся ближе всего к началу координат. действий.
Неформально следящий алгоритм может быть описан следующей последовательностью 1. Первоначально выбирается следящая кубическая область с длиной ребра r и нулевой вершиной в точке q, заведомо покрывающая многогранник. 2. Выбирается целевая точка z = T c, расположенная вне следящей области. 3. Следящая область разбивается на K ячеек. 4. В условиях динамического изменения входных данных (A, b, c), для всех ячеек внутри следящей области вычисляется псевдопроекция из точки z на пересечение ячейки с многогранником. Если пересечение пусто, то такие ячейки отбрасываются. щей области в w раз и переходим на шаг 2. 6. Если получено непустое множество псевдопроекций, то на нем вычисляется максимум целевой функции.
ребра следящей области r уменьшается в 2 раза. 7. Если расстояние от точки максимума до центра следящей области меньше 41 r, то длина ребра следящей области r увеличивается в 2 раза. 8. Если расстояние от точки максимума до центра следящей области больше 34 r, то длина 9. Центр следящей области перемещается в центр ячейки, в которой достигнут максимум, и переходим на шаг 2.
Псевдопроекции на шаге 4 для различных ячеек следящей области могут вычисляться параллельно без обменов данными между MPI-процессами. Фейеровский процесс вычисления псевдопроекции, стартовавший из целевой точки на шаге 4, вычисляется до конца, несмотря на то, что координаты целевой точки и сам многогранник в процессе счета могли измениться. Выполнение условия шага 5 означает, что входные данные изменяются быстрее, чем мы вычисляем псевдопроекции. Константы 1/2 и 3/4 , используемые при выполнении шагов 7 и 8, являются значениями параметров алгоритма. 4. Межузловое распараллеливание
Межузловое распараллеливание – это распараллеливание вычислений между отдельными узлами кластерной вычислительной системы. В качестве технологии параллельного программирования используется MPI. Каждый MPI-процесс работает на отдельном узле кластера. Обмен данными между MPI-процессами осуществляется с помощью передачи сообщений по коммуникационной сети. Один MPI-процесс считает одну псевдопроекцию на пересечение одной ячейки следящей области с многогранником. Таким образом, количество ячеек следящей области всегда равно константе K, совпадающей с количеством MPIпроцессов.
Опишем метод межузлового распараллеливания следящего алгоритма. Без ограничения общности мы можем считать, что все процессы происходят в положительной области координат. Пусть P – количество доступных MPI-процессов (количество процессорных узлов в многопроцессорной системе). Определим общее количество ячеек в следящей области следующим образом:
Тогда количество ячеек, примыкающих к одному ребру куба следящей области, равно h = K1/n. Зададим в пространстве целочисленных координат u0, . . . , un линейную нумерацию ячеек следящей области следующим образом. Пусть ячейка \alpha имеет целочисленные координаты (\alpha 0, . . . , \alpha n). Тогда ее номер k\alpha вычисляется по формуле:
K =
B\igl(for
P 1/n\Bigrflo) n
. k\alpha = \sum Рис. 1. Линейная нумерация ячеек следящей области при n = 3. На рис. 1 приведен пример такой линейной нумерации при n = 3. Например, ячейка с номером 19 на рис. 1 имеет целочисленные координаты (1, 0, 2). Действительно, 19 = 1 \cdot 30 + 0 c\dot 31 + 2 c\dot 32.
Выразим целочисленные координаты ячейки \alpha через ее порядковый номер k\alpha . Из (
\alpha 0 = k\alpha mod n; \alpha 1 = k\alpha - \alpha 0 mod n;
n \alpha 2 = k\alpha - \alphan02- \alpha 1n mod n;
. . . . . . . . . \alpha i =
i- 1 k\alpha - \sum \alpha jnj j=0 ni
mod n.
\alpha 0 = k\alpha - (k\alpha d\iv n) \cdot n.
Таким образом, в общем виде имеем:
Формула (9) содержит ресурсоемкую операцию возведения в степень. От нее можно
избавиться следующим образом. По определению операции mod из (
\alpha 1 = k\alpha d\iv n - ((k\alpha d\iv n) \div n) \cdot n. Подставив в (8) вместо \alpha 0 правую часть уравнения (10), а вместо \alpha 1 – правую часть уравнения (12), получим \alpha 2 = k\alpha - \alphan02- \alpha 1n mod n = k\alpha - (k\alpha - (k\alpha d\iv n)\cdot n)- (nk\2alpha d\iv n- ((k\alpha d\iv n)\div n)\cdot n)c\dot n mod n = ((k\alpha d\iv n) \div n) mod n. Из (11) и (13) для i = 1, . . . , n - 1 по индукции получаем: \alpha i = (((k\alpha d\iv n) . . .) \div n) modn.
i
u\nderbac{} Пусть g = (g0, . . . , gn- 1) – нулевая вершина куба следящей области; y = (y0, . . . , yn- 1) – нулевая вершина произвольной ячейки \alpha . Выразим координаты точки y через координаты r точки g. Обозначим s = K1/n – шаг сетки. Тогда для i = 0, . . . , n - 1. тами (\gama 0, . . . , \gama n- 1), где
Определим в качестве центральной ячейки куба ячейку \gama с целочисленными координаточки q через координаты точки g, используя формулу (15): Пусть q = (q0, . . . , qn- 1) – нулевая вершина центральной ячейки \gama . Выразим координаты u\nderbac{} l\eft{
yi = gi + s\alpha i \gama 0 = . . . = \gama n- 1 =
B\iglfor
K1/n\Big/ 2
B\igrflo qi = gi + s\gama i, - x0 \leq c\dot \cdot \cdot c\dot c\dot \cdot c\dot \cdot c\dot c\dot \cdot c\dot
y0 - xn- 1 \leq
yn- 1 x0 \leq y0 + s c\dot \cdot \cdot c\dot c\dot \cdot c\dot \cdot c\dot
c\dot \cdot c\dot xn- 1 \leq yn- 1 + s для i = 0, . . . , n - 1. задается системой из 2n неравенств:
Пусть y – нулевая вершина ячейки \alpha .Тогда область внутри ячейки \alpha (включая границы) неравенств с n неизвестными. Чтобы найти пересечение ячейки \alpha с многогранником M необходимо к системе неравенств Ax \leq b добавить неравенства из системы (18). При этом получается система из m + 2n Параллельная версия следящего алгоритма изображена на рис. 2 в виде диаграммы деятельности UML. В качестве начального значения длины ребра r кубической следящей области берется значение R, а в качестве нулевой вершины g кубической следящей области берется точка G, которые заведомо обеспечивают покрытие многогранника M следящей областью. Начальный шаг сетки s определяется по формуле s = K1r/n , где количество ячеек K следящей области вычисляется по формуле (4). Координаты целевой точки z получаются путем умножения вектора c коэффициентов целевой функции на масштабирующий коэффициент T , выбираемый таким образом, чтобы целевая точка находилась вне следящей области. Целочисленные координаты центральной ячейки следящего куба определяется по формуле (16). Нулевая вершина q центральной ячейки следящей области соответственно вычисляется по формуле (17).
Вектор-функция \pi (z, k) с помощью фейеровского процесса, описанного в разделе 2, вычисляет псевдопроекцию целевой точки z на пересечение многогранника с ячейкой, имеющей порядковый номер k, и возвращает точку псевдопроекции xk или вектор (- 1, . . .), (14) (15) (16) (17) (18) если точка псевдопроекции не принадлежит многограннику (пересечение многогранника M с ячейкой пусто).
Псевдопроекции вычисляются параллельно без обменов данными в цикле until. В ходе вычисления псевдопроекций исходные данные задачи (1) могут меняться. После того как все точки псевдопроекций x0, . . . , xK- 1 вычислены, в цикле for происходит вычисление порядкового номера k\alpha ячейки \alpha , на которой достигается максимум целевой функции в точке псевдопроекции. При этом в переменной C вычисляется максимум целевой функции. В соответствии с этим в качестве начального значения k\alpha выбирается наименьшее целое число в машинном представлении MinInt, а в качестве начального значения C выбирается наименьшее вещественное число в машинном представлении MinFloat.
После завершения цикла for проверяется значение k\alpha . Если оно равно MinInt, значит пересечение всех ячеек с многогранником M оказалось пустым. Это означает, что в результате изменения исходных данных задачи (1) многогранник M «уплыл» за пределы следящей области. В этом случае в w раз увеличиваются: длина r ребра следящей области, шаг сетки s и координаты целевой точки z. Константа w является параметром алгоритма.
Если значение k\alpha отлично от MinInt, значит пересечение многогранника M со следящей областью не пусто. В этом случае в качестве новой центральной ячейки рассматривается ячейка с номером k\alpha и вычисляются координаты q\prime – нулевой вершины этой ячейки. Для этого используется вектор-функция zero(k\alpha ), диаграмма деятельности которой приведена на рис. 3.
Рис. 3. Вектор-функция, вычисляющая нулевую вершину ячейки с порядковым номером k\alph .
Далее рассматриваются три случая в зависимости от удаленности нулевой вершины новой центральной ячейки следящей области от старой. Если \| q - q\prime \| > 43 r, то r, s и координаты z увеличиваются в полтора раза. Если \| q - q\prime \| < 41 r, то r, s и координаты z уменьшаются в полтора раза. Если 14 r \leq \| q - q\prime \| l\eq 43 r то r, s и координаты z остаются прежними. После этого происходит сдвиг нулевой вершины следящей области на вектор (q\prime - q) и копирование координат q\prime в q. Далее выполняется очередная итерация цикла until. 5. Заключение
В работе была описана параллельная реализация следящего алгоритм для решения нестационарных задач большой размерности на кластерных вычислительных системах. Алгоритм был реализован на языке Си с использованием библиотеки MPI. Начальные эксперименты на модельной задаче показали почти линейную масштабируемость параллельной реализации следящего алгоритма на 480 узлах суперкомпьютера «Торнадо ЮУрГУ». Авторы планируют реализовать внутриузловое распараллеливание следящего алгоритма с использованием технологии параллельного программирования OpenMP. Объектом распараллеливания на этом уровне является отдельный фейеровский процесс, вычисляющий псевдопроекцию. Для каждого подвектора организуется отдельная нить. Вычисления предполагается выполнять с использованием ускорителей Xeon Phi, которыми оснащен вычислительный кластер «Торнадо ЮУрГУ». Эффективное использование многоядерных ускорителей достигается путем распараллеливания процесса вычисления отдельной псевдопроекции. При этом предполагается использовать метод разбиения вектора на подвекторы [9]. Суть метода заключается в том, что вектор исходной точки делится на подвекторы по числу процессорных ядер многоядерного ускорителя. На каждом подвекторе независимо делается v фейеровских приближений с использованием редуцированных фейеровских отображений. Такое редуцированное отображение воздействует только на «свой» подвектор, оставляя оставшуюся часть вектора без изменений. Затем нити управления через общую память обмениваются модифицированными подвекторами и процесс продолжается. В работе [9] была доказана сходимость описанного итерационного процесса. Литература Solving unstable linear programming problems of high dimension on cluster computing systems Irina Sokolinskaya and Leonid Sokolinsky Keywords: linear programming, unstable problems of high dimension, cluster computing systems, Fejer mappings