Одной из особенностей современного экономико-математического моделирования является появление нестационарных задач линейного программирования [1] большой размерности с быстро меняющимися входными данными. Кроме этого, часть ограничений может оказаться плохо формализуемой [2]. Одним из примеров таких задач является задача управления портфелем ценных бумаг с использованием методов алгоритмической торговли [3,4]. В подобных задачах количество переменных и неравенств в системе ограничений может составлять десятки и даже сотни тысяч, а период изменения исходных данных находиться в пределах сотых долей секунды. В работе [5] был описан параллельный алгоритм для решения задач линейного программирования с плохо формализуемыми ограничениями. Для преодоления проблемы нестационарности входных данных в работах [6,7] был предложен «следящий» алгоритм решения задачи линейного программирования с использованием фейеровских отображений [8], ориентированный на кластерные вычислительные системы с многоядерными ускорителями. В данной работе дается параллельная реализация следящего алгоритма с использованием диаграммы деятельности UML. Статья организована следующим образом. В разделе 2 приводится формальная постановка задачи линейного программирования, даются определения фейеровского процесса и операции псевдопроектирования на многогранник. В разделе 3 приводится неформальное описание следящего алгоритма. В разделе 4 описывается межузловое распараллеливание следящего алгоритма, предполагающее использование технологии параллельного программирования MPI. В заключении суммируются полученные результаты и определяются направления дальнейших исследований.

Пусть задана задача линейного программирования max \{ \langle c, x\rangle | Ax \leq b, x \geq 0\} .

(

Определим отображение \varphi : \BbbR n \rightar \BbbR n следующим образом:

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Пусть M -многогранник, задаваемый ограничениями задачи линейного программирования (1). Такой многогранник всегда является выпуклым. Известно [8], что \varphi будет однозначным непрерывным M -фейеровским отображением для любой системы положительных коэффициентов \{ \alpha j > 0\} (j = 1, . . . , m), таких, что m \sum j=1 \alpha j = 1, и коэффициентов релаксации 0 < \lambda j < 2. Полагая в формуле (2) \lambda j = \lambda и \alpha j = 1/m (j = 1, . . . , m), получаем формулу

В общем виде следящий алгоритм может быть описан следующим образом. Все пространство \BbbR n делится гиперкубической сеткой на ячейки с переменной длиной грани. Размер ячейки может динамически меняться в зависимости от скорости изменения исходных данных: чем быстрее меняются данные, тем больше становится ячейка, чтобы успевать следить. Алгоритм постоянно находит (отслеживает) ячейку наименьшего размера, на которой достигается максимум целевой функции. Будем такую ячейку называть целевой. В каждой итерации алгоритма просчитываются ячейки, входящие в некоторую следящую область вокруг целевой ячейки. В простейшем случае следящая область может быть кубической формы со следящей ячейкой в центре. В общем случае следящая область может быть произвольной выпуклой фигурой. При этом ячейки всегда имеют кубическую форму. Целевая точка -это некоторая точка, находящаяся вне следящей области. Ее координаты могут динамически вычисляться по коэффициентам целевой функции. Например, при целевой функции Q (x) = x 1 + 2x 2 в качестве целевой точки можно взять точку (T, 2T ), где Tдостаточно большое положительное число. Нулевой вершиной кубической области будем называть вершину куба, находящуюся ближе всего к началу координат. Нулевой вершиной кубической ячейки будем называть вершину ячейки, находящуюся ближе всего к началу координат.

Неформально следящий алгоритм может быть описан следующей последовательностью действий.

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

2. Выбирается целевая точка z = T c, расположенная вне следящей области.

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

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

Тогда количество ячеек, примыкающих к одному ребру куба следящей области, равно h = K 1/n . Зададим в пространстве целочисленных координат u 0 , . . . , u n линейную нумерацию ячеек следящей области следующим образом. Пусть ячейка \alpha имеет целочисленные координаты (\alpha 0 , . . . , \alpha n ). Тогда ее номер k \alpha вычисляется по формуле:

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org

. . . . . . . . .

Формула ( 9) содержит ресурсоемкую операцию возведения в степень. От нее можно избавиться следующим образом. По определению операции mod из (6) получаем

С помощью символа \div здесь обозначается целочисленное деление. Подставив в (7) вместо \alpha 0 правую часть этого уравнения, получим

По определению операции mod отсюда следует

Подставив в (8) вместо \alpha 0 правую часть уравнения (10), а вместо \alpha 1 -правую часть уравнения (12), получим

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Из (11) и (13) для i = 1, . . . , n -1 по индукции получаем:

Пусть g = (g 0 , . . . , g n - 1 ) -нулевая вершина куба следящей области; y = (y 0 , . . . , y n - 1 )нулевая вершина произвольной ячейки \alpha . Выразим координаты точки y через координаты точки g. Обозначим s = r K 1/n -шаг сетки. Тогда

Определим в качестве центральной ячейки куба ячейку \gamma с целочисленными координатами (\gamma 0 , . . . , \gamma n - 1 ), где

Пусть q = (q 0 , . . . , q n - 1 ) -нулевая вершина центральной ячейки \gamma . Выразим координаты точки q через координаты точки g, используя формулу (15):

Пусть y -нулевая вершина ячейки \alpha .Тогда область внутри ячейки \alpha (включая границы) задается системой из 2n неравенств: \left\{ - x 0 \leq - y 0 \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot - x n - 1 \leq - y n - 1

x 0 \leq y 0 + s \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot

Чтобы найти пересечение ячейки \alpha с многогранником M необходимо к системе неравенств Ax \leq b добавить неравенства из системы (18). При этом получается система из m + 2n неравенств с n неизвестными. Параллельная версия следящего алгоритма изображена на рис. 2 в виде диаграммы деятельности UML. В качестве начального значения длины ребра r кубической следящей области берется значение R, а в качестве нулевой вершины g кубической следящей области берется точка G, которые заведомо обеспечивают покрытие многогранника M следящей областью. Начальный шаг сетки s определяется по формуле s = r K 1/n , где количество ячеек K следящей области вычисляется по формуле (4). Координаты целевой точки z получаются путем умножения вектора c коэффициентов целевой функции на масштабирующий коэффициент T , выбираемый таким образом, чтобы целевая точка находилась вне следящей области. Целочисленные координаты центральной ячейки следящего куба определяется по формуле (16). Нулевая вершина q центральной ячейки следящей области соответственно вычисляется по формуле (17).

Вектор-функция \pi (z, k) с помощью фейеровского процесса, описанного в разделе 2, вычисляет псевдопроекцию целевой точки z на пересечение многогранника с ячейкой, имеющей порядковый номер k, и возвращает точку псевдопроекции x k или вектор ( - 1, . . .), если точка псевдопроекции не принадлежит многограннику (пересечение многогранника M с ячейкой пусто).

Псевдопроекции вычисляются параллельно без обменов данными в цикле until. В ходе вычисления псевдопроекций исходные данные задачи (1) могут меняться. После того как все точки псевдопроекций x 0 , . . . , x K - 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 \alpha .

Далее рассматриваются три случая в зависимости от удаленности нулевой вершины новой центральной ячейки следящей области от старой. Если \| q -q \prime \| > 3 4 r, то r, s и координаты z увеличиваются в полтора раза. Если \| q -q \prime \| < 1 4 r, то r, s и координаты z уменьшаются в полтора раза. Если 1 4 r \leq \| q -q \prime \| \leq 3 4 r то r, s и координаты z остаются прежними. После этого происходит сдвиг нулевой вершины следящей области на вектор (q \prime -q) и копирование координат q \prime в q. Далее выполняется очередная итерация цикла until.

В работе была описана параллельная реализация следящего алгоритм для решения нестационарных задач большой размерности на кластерных вычислительных системах. Алгоритм был реализован на языке Си с использованием библиотеки MPI. Начальные эксперименты на модельной задаче показали почти линейную масштабируемость параллельной реализации следящего алгоритма на 480 узлах суперкомпьютера «Торнадо ЮУрГУ». Авторы планируют реализовать внутриузловое распараллеливание следящего алгоритма с использованием технологии параллельного программирования OpenMP. Объектом распараллеливания на этом уровне является отдельный фейеровский процесс, вычисляющий псевдопроекцию. Для каждого подвектора организуется отдельная нить. Вычисления предполагается выполнять с использованием ускорителей Xeon Phi, которыми оснащен вычислительный кластер «Торнадо ЮУрГУ». Эффективное использование многоядерных ускорителей достигается путем распараллеливания процесса вычисления отдельной псевдопроекции. При этом предполагается использовать метод разбиения вектора на подвекторы [9]. Суть метода заключается в том, что вектор исходной точки делится на подвекторы по числу процессорных ядер многоядерного ускорителя. На каждом подвекторе независимо делается v фейеровских приближений с использованием редуцированных фейеровских отображений. Такое редуцированное отображение воздействует только на «свой» подвектор, оставляя оставшуюся часть вектора без изменений. Затем нити управления через общую память обмениваются модифицированными подвекторами и процесс продолжается. В работе [9] была доказана сходимость описанного итерационного процесса. Литература