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

УДК 519.852

ПАРАЛЛЕЛЬНЫЕ РЕАЛИЗАЦИИ СИМПЛЕКС-МЕТОДА

ДЛЯ БЕЗОШИБОЧНОГО РЕШЕНИЯ

ЗАДАЧ ЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ

А.В. Панюков, В.В. Горбик

THE PARALLEL SIMPLEX-METHOD ACHIEVEMENTS

FOR ERRORLESS SOLVING

OF LINEAR PROGRAMMING PROBLEMS

A.V. Panyukov, V.V. Gorbik

В работе рассмотрены подходы к решению задачи линейного программирования с абсолютной точностью, достигаемой применением в алгоритмах симплекс-метода дробно-рациональных вычислений без округления. Если при этом т - минимальная из размерностей задачи, I - число бит, необходимых под один численный элемент исходных данных, то пространственная сложность алгоритма не превосходит 4/т4 + о(т3), при этом вычислительная сложность одной итерации симплекс-метода не превосходит 0(/т4), а эффективность распараллеливания (т.е. отношение ускорения к числу процессоров) в предложенной реализации параллельного алгоритма составляет в асимптотике 100%.

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

Techniques of obtaining exact solutions of linear programming problems are subjects of this paper. Absolute accuracy are arrived at implementation of simplex-algorithm with exact rational-fractional computation. In this case if m is minimal of problem dimensions, and / is number of bits for a source data item then space complexity are no more 4/m4 + o(rn3), one iteration time complexity are no more 0(lrn4), and paralleling efficiency (i.e. ratio of acceleration to number of processors) asymptotical estimate are 100%.

Keywords: linear programming, simplex method, distributed computing, parallel computing, rational computations, optimization, arbitrary precision, interval arithmetic.

Введение

В настоящее время широко распространены предрассудки, не основанные на доказательствах и порождающие ошибки в расчетах: (1) распространение свойства ассоциативности операций сложения и умножения в поле действительных чисел на конечное множество машинных «действительных> чисел; (2) распространение свойства непрерывной зависимости от параметров решения системы, полученные после «:эквивалентных> преобразований, на исходную систему имеют популярные коммерческие пакеты MatLab, MathCad и т.п., а также свободно распространяемый пакет SciLab. Использование в вычислениях разного числа процессоров во многих случаях дает существенно различающиеся результаты, демонстрируя необходимость доказательных вычислений (см. работы [1, 2]).

Потенциал имеющихся пакетов, поддерживающих символические вычисления, не позволяет решать реальные проблемы математического и имитационного моделирования. Возможность обеспечения вычислений с произвольной точностью в программах пользователя дает библиотека GMP. Однако, ее использование требует от пользователя разработки собственного интерфейса для организации распределенных и параллельных вычислений [3]. Развитием библиотеки GMP является библиотека ExactComputational [4], которая предоставляет своим объектам возможность их использования в параллельных вычислениях.

Анализ способов распараллеливания показывает эффективность распараллеливания «по информации:». При этой технологии вычислительный процесс строится на основе единственной программы, запускаемой на всех процессорах вычислительной системы или на многих станциях локальной сети. Копии программы могут выполняться по разным ветвям алгоритма, обрабатывая подмножества данных. Неизбежна синхронизация во времени и при обработке общих данных. Данная идеология используется в стандарте MPI (Message Passing Interface ).

Ранее авторами разработаны алгоритмы и программное обеспечение для абсолютно точного решения систем линейных алгебраических уравнений [5] и вычисления обобщенной обратной матрицы Мура-Пенроуза [6] на многопроцессорных вычислительных системах с использованием классов overlong и rational из библиотеки ExactComputational [4]. Теоретическое и практическое исследование данного программного обеспечения демонстрирует высокую эффективность использования многопроцессорных вычислительных систем.

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

1. Техника реализации симплекс-метода

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

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

max {стх : Ах = b > 0, х > 0; с, х £ Rn; b £ Rm} . (1)

1.1. Метод симплекс-таблиц

Данный метод на итерации к пересчитывает симплекс-таблицу

z(k) = ^В(к)в^)~1ъ z(k) = -ст + е^Б^) гА

Хв(к) = В(к)-гЬ В{к)~1А

где В (к) - базисная матрица, содержащая все относящиеся к базисным переменным к-й итерации столбцы матрицы А (базисные столбцы); св^ - вектор коэффициентов целевой функции, относящихся к базисным переменным к-й итерации. При этом левый столбец симплекс-таблицы содержит вектор X= В(к)~1Ь значений базисных переменных к-й итерации и значение г ^ целевой функции на этом решении. Столбцы матрицы 5, соответствующие базисным переменным, являются ортами (т.е. из них может быть составлена

единичная матрица). Верхняя строка содержит вектор = —с? + с^^В(к) 1А невязок

двойственной задачи. Значения элементов строки г(к), относящиеся к базисным столбцам являются нулевыми.

Критерием оптимальности текущего базисного решения является неотрицательность

(к)

вектора г. Если же существует небазисная переменная х{ : г\ < 0, то условие Ь =

{I : 3(к) > 0} = 0 является признаком неограниченности целевой функции. Если же

Ь = {I : 5(к) > 0} ф 0, то введение переменной х* в число базисных приведет к увели-

чению целевой функции на величину

в(*0г*

(к)

С

(к)

1*1

где I* = аг§1ехтіп

ІЄЬ С

(к)

и

(2)

Данный способ выбора I* сохранит неотрицательность базисного решения на следующей итерации и исключает зацикливание, когда А* = 0 [7].

Переход от таблицы к-й итерации к таблице (к + 1)-й итерации осуществляется применением процедуры исключения Жордана-Гаусса для столбца г (ведущего столбца), используя в качестве ведущей строку I*

с(к+1) ___ ,

о(к)

С

С

С

(к) і* з

(к)

1*3 о(к) (к) ,

1*г

о(к) ’ й1*і

если І ф I*,

если 1 = 1*

(3)

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

1.2. Метод обратной матрицы

В данном методе, в отличие от табличного, на каждой итерации вместо пересчета симплекс-таблицы пересчитывается матрица В (к)-1 . Наличие обратной матрицы позволяет для текущего базиса легко находить у{к)Т = с^^іЗ(£;)-1 - соответствующее решение двойственной задачи. Текущий базис является оптимальным если соответствующее ему двойственное решение допустимо: у(к)тА > ст. Если же в матрице А найдется столбец Аі(к) : у(к)тАі(к) < с,£(к), то введение его в число базисных приведет к увеличению целевой функции.

Образом столбца в симплекс-таблице 3(к) будет вектор д = Вік)-1 А^уу Поэтому ведущей строкой, определяющей выводимый из базиса столбец, будет

X

г = аг§1ех тіп

"в(к)і

І: уі>0

а новые значения базисных переменных равны

/ Г х

9і

(УІ = 1, 2,..., ш)

В(к+1)і \ Хв{к)і

91 у

в(к)і ~ -АВДг1

9г =

если І Ф г, ^ если І = г

(4)

(5)

Базисные матрицы В (к) = иВ(^+1) = ( ) отличаются только г-м столб-

V э /1^=1 V э /г,5=1

цом, поэтому элементы обратной матрицы В (к + I)-1 = \$и+1) следующим образом

V ■' /г, 1=1

(6)

Оценим алгебраическую вычислительную сложность рассматриваемого метода. Очевидно, что в общем случае на заключительной итерации потребуется проверка всех п ограничений двойственной задачи, для чего потребуется не более тп операций умножения и т(п — 1) операций сложения. На промежуточных итерациях для установления недонусти-мости двойственного решения эта величина, как правило, является существенно меньшей. Пересчет обратной матрицы потребует не более т операций деления, не более т2 операций умножения и не более (т2 — 1) операций сложения/вычитания. Поскольку т < п,то затраты вычислительных ресурсов на пересчет обратной матрицы могут быть существенно ниже затрат на пересчет симплекс-таблицы. Алгебраическая пространственная сложность метода обратной матрицы определяется числом элементов в исходных данных и в обратной матрице, т.е. равна тп + т2 + О(т), т.е. больше чем в табличном симплекс-методе.

1.3. Битовая сложность абсолютно точной реализации симплекс-метода

Практическая реализуемость вычислений без округления, в частности, безошибочного решения задачи линейного программирования, определяется требуемыми для вычислений ресурсами: числом бит оперативной памяти и количеством операций с битами. Далее через 5'(Л) будем обозначать число бит, требуемое для представления объекта Л, а через С(Х) - число битовых операций, выполненных при нахождении представления объекта Л. Например, число бит, требуемое для представления целых чисел, будет равно 5(0) = 1, (Уг 6 Z \ {0})5'(г) = |"1о§2 |г|]. Число бит, требуемое для представления рационального числа г = р)д, р, ^/eZ имеет верхнюю оценку Б (г) < 0(Б(р) + Я(д)).

Легко проверить, что если Б(р), Б(д) - память, требуемая для представления рациональных чисел р, д, то память, требуемая для представления результата арифметической операции о £ {+, —,/, х} над данными числами будет Б(р од) < 5'{р) + 8{д). Для битовой вычислительной сложности выполнения операции о £ { + , —, /, х} с помощью классических арифметических алгоритмов (умножение/деление столбиком) справедлива оценка С(р о д) < О (5'(р)5'(д)). Использование алгоритмов быстрого умножения дает оценку С(р° я) < (&(р) + 3(д))> которая будет использована в работе.

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

Предложение 1. Пусть В = (Ьц) - целочисленная гпхгп матрица, I = тах 5'(Ь^)-

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

і^=1,2,...,п

Тогда

В) < п(\о^2 т + I)

т

^ В\ = ^2 П Ька(к) < ^ П \Ька(к) I < т\Ьгп < (тЬ)т

о- к=1

с к=1

где L = max{|bjj | : г, j = 1, 2,..., m}.

Из данной оценки следует S'(det В) = log2 | det В\ < rn(log2 гп + 1). □

Предложение 2. Пусть В = (bij) - т х т рациональная матрица, I = max S(bij). Тогда S'(det В) < m(log2 m + (2m + 1)/).

Доказательство. Пусть Kr равно наименьшему общему кратному всех знаменателей строки г = 1,2,... ,т матрицы В. Очевидно, что S(Kr) < 1т. Пусть К представляет диагональную матрицу: diag(K) = {Kr : г = 1, 2,...,m}. Рассмотрим матрицу В = КБ. Данная матрица является целочисленной. Верхняя оценка числа бит, необходимых для одного элемента матрицы В, имеет вид

I = max S(bij) = max S(Kibij) < l{rn + 1). г,2=1,2,...,m г,2=1,2,...,m

Принимая во внимание утверждение 1, а также log2 | det К~1 \ < 1т, получим S'(det В) = S ^det К-1 ■ det(i3)^ < m(log2 rn + (2m + 1)/)

□

Из формул Крамера для систем линейных алгебраических уравнений и доказанных утверждений вытекает

Следствие 1. Если один численный элемент исходных данных имеет пространственную сложность не более I, то

• элементы симплекс-таблицы и обратной матрицы имеют пространственную сложность не более 4/m2 + Zm + m log2 m + 1;

• столбец симплекс-таблицы и обратной матрицы имеют пространственную сложность не более 4Zm3 + Zm2 + m2 log2 m + m;

• для представления симплекс-таблицы достаточно Alm^n + о(1т2п) бит;

• для представления обратной матрицы достаточно 41т4 + o(Zm3) бит.

Поскольку m < п, то, относительно битовой пространственной сложности симплекс-метода, предпочтение следует отдать методу обратной матрицы.

Для анализа вычислительной сложности и эффективности распараллеливания симплекс-метода рассмотрим его параллельные версии.

2. Параллельные версии симплекс-метода

2.1. Метод симплекс-таблиц

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

Принятые соглашения позволяют предложить следующую параллельную реализацию итерации метода симплекс-таблиц.

Процесс К = 1, 2,..., N

Боо = % г\(К-ы1)п^+2 • • • Хг Кп і 1 ^ 1

Біо = Хв і Б20 = ХВ 2 Бт0 ХВт ^т\(кк1)пї+1 3і\(Км1)пї+2 б2|-(*--і)гі-|+2 *^тГ(К «1)п1+2 • • • Біг Кп і 1 N 1 Б2Г Кп і 1 N 1 БтГ Кпп 1 N 1

Рис. 1. Декомпозиция симплекс-таблицы по процессорам

Алгоритм ТаЬиІагБітрІех Итерация к

• Данные: симплекс-таблица в(К) для каждого процесса К = 1,2,3,...,

• Шаг 1. Каждому процессу К = 1,2,3,..., -/V

— найти столбец ік : г{к < 0;

— если столбец і к не найден, то положить С (К) = А (К) = і к = 0 и перейти на шаг 2;

— найти строку

ХвЛ Энк)'

— если строка 1к не найдена, то завершить выполнение всех процессов и вернуть ”Не ограничена”;

— вычислить Аік = —ХвікХік/вікік1 положить А(К) = Аік, С(К) = %к и перейти на шаг 2.

Комментарий. При выполнении данного шага будет либо установлена неразрешимость задачи, либо найдены данные для изменения базиса: ведущий столбец ік -кандидат для ввода в базис, ведущая строка І к, определяющая столбец выводимый из базиса, и приращение целевой функции А{к, - либо установлено отсутствие таких кандидатов: і к = 0.

• Шаг 2. Для Ь = 1, 2, 4, 8,..., N каждому процессу К = 1, 2, 3,..., N, номер которого удовлетворяет условию {{К — 1)%(2Ь)) < Ь, осуществить обмен данными с процессом К + Ь:

— Если С (К) = 0, то положить С (К) = С (К + Ь), А (К) = А (К + Ь) и продолжить вычисления для следующего Ь.

— Если С (К + V) = 0, то положить С (К + V) = С (К), А (К + Ь) = А (К) и продолжить вычисления для следующего Ь.

— Если А (К) > А (К + Ь), то положить С (К + Ь) = С (К), А (К + Ь) = А (К), иначе положить С (К) = С {К + I), А (К) = А (К + I).

— Продолжить вычисления для следующего Ь.

Ік = тіп

І' $ІІК >0

Комментарий. После завершения данного шага в каждом процессе К = 1, 2, 3,..., N значение С (К) будет равно номеру ведущего процесса

К* =

аг§1ех тах А*к, если {К : гк Ф &} Ф К\ гкфО

О,

в противном случае,

Шаг З.Если К* = 0, то каждому процессу К = 1, 2, 3,..., N для

Г(К — 1)п_. Г(К — 1)и, г Кп..

*=гЧИ-1 + 1. 1

ПОЛОЖИТЬ

Хвь если («З'а- = 1) Д {(уг = о, 1,... ,1 - 1,1 + 1,... ,т) в,1к = 0),

О, в противном случае.

Вернуть х - оптимальное решение задачи, Е - оптимальное значение целевой функции и завершить выполнение алгоритма.

Комментарий. Если текущее базисное решение задачи оптимально, то формирется ответ и завершается выполнение алгоритма.

• Шаг 4. Ведущему процессу К* передать остальным процессам ведущий столбец

\Т

и номер ведущей строки 1к* ■

• Шаг 5. Каждому процессу К = 1, 2, 3,..., N пересчитать по формуле (3) симплекс-таблицу £(К). Перейти к следующей итерации.

• Конец алгоритма

Из изложенного выше видно, что описание параллельной версии табличного симплекс-метода не намного сложнее чем для обычного алгоритма. При этом используется 1о§2 ^] широковещательных коммуникаций между процессами при нахождении ведущего процесса К* и одна широковещательная коммуникация для передачи ведущего столбца 5^* и числа 1к* ■ В таблице 1 приведена сводка требуемого алгебраического (т.е. количества используемых алгебраических операций) вычислительного ресурса для последовательной и параллельной реализаций метода симплекс-таблиц. Из табл. и описания алгоритма видно, что процессоры загружены равномерно, а эффективность распараллеливания растет с ростом сложности задачи, достигая в пределе 100%.

2.2. Метод обратной матрицы

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

Л(К ) =

С(К)Т = ( С|-Й^1+1 сГЙ£_^1+2 ( +1) а2(г^1+1) а2(г^1+2)

С\Кп] ), К = 1,2,

04(^1) \

а2(Г^)

N;

■1+1

ат(,{к_12~1+2\

1№) )

(7)

К = 1, 2, 3, ...^. (8)

а

N

Таблица 1

Алгебраический вычислительный ресурс реализаций метода симплекс-таблиц

Один процессор N процессоров

Количество Количество Нагрузка

Оператор алгебраических пересылаемых на один

операций операндов процесс

Проверка условия оптимальности [1,п] - [1,п^ ]

Определение ведущей строки 2т + п + 2 - 2т + n/N + 2

Выбор ведущего процесса - N Г1оё2 N1

Пересылка ведущего столбца - т + 2 -

Пересчет симплекс-таблицы 2(т + 1)(п + 1) - 2т(1 + n/N)

Итого: [1, п] + 2тп+ 4(т + п +1) т + N + 2 [1, п/^ ] + 2(n/N )(т + 1) + 4т + 2 + 1о§2 N

При этом предполагается, что вектор двойственных переменных у размещен в каждом процессе К = 1, 2,..., N,

Эффект от распараллеливания при пересчете обратной матрицы В-1 будет максимальным при ее декомпозиции по строкам в равных пропорциях на блоки по числу процессов

( Ъ

в~\ку

{К- 1)71

1+і) 1 Ь(Г^^1+1)2 Ь(г^^1+і)т ^

• Ь(г(к^ї1+2)т

Ъ^(к-1)т^+2у Ь^(к-1)гг^+2у

V ь(г^^^1)1

ъ

(г ^1 +2)2

ъ

К = 1, 2, З, ...,^ (9)

(Г^Р1+2)т /

Из описания метода обратной матрицы следует целесообразность декомпозиции вектора Хв базисных переменных, вектора Св коэффициентов целевой функции при базисных переменных и вектора д по таким же блокам строк

Хв (К) =

( Хв\(к-м1)т^+і \ Хв1(К^Г)гп^+2

V Хв г 3р1

; св (К) =

( Свг(^1+1 ^ Св\(К^Г)гп^+2

V Свг^1 /

; д(к) =

5гЙ£_і)^і+2

V /

К = 1, 2, З, ...,^ (10)

Способ (9) декомпозиции обратной матрицы по процессам позволяет каждому процессу вычислить субвектор двойственных переменных

у(К) = (В(К)~1')Т сВ(К)-

Очевидно, что вектор двойственных переменных

(11)

N

К=1

Изложенное выше позволяет предложить следующую параллельную реализацию итерации метода обратной матрицы.

N

Алгоритм 1^ВазМа1г81тр1ех Итерация к

• Данные: в каждом процессе К = 1, 2, 3,..., N размещены

у - вектор двойственных переменных, построенный по текущему базисному решению прямой задачи;

М(К) - вектор номеров базисных переменных процесса;

Хв{К) - вектор значений базисных переменных процесса;

с в (К) - вектор значений коэффициентов целевой функции при базисных переменных процесса;

А(К) - блок матрицы ограничений процесса;

с(К)т - блок коэффициентов целевой функции процесса;

В-1 (К) - блок обратной матрицы процесса.

• Шаг 1. Каждому процессу К = 1,2,3,..., -/V

если столбец г к не найден, то положить С (К) = А (К) = гк = 0; иначе положить

• Шаг 2. Для Ь = 1, 2,4, 8,..., N каждому процессу К = 1,2,3,..., N, номер которого удовлетворяет условию ((К — 1)%(2Ь)) < Ь, осуществить обмен данными с процессом

Если С (К) = 0, то положить С (К) = С (К + Ь), А (К) = А (К + Ь) и продолжить вычисления для следующего Ь.

Если С (К + Ь) = 0, то положить С(К + Ь) = С (К), А (К + V) = А (К) и продолжить вычисления для следующего Ь.

Если А (К) < А (К + Ь), то положить С (К + Ь) = С (К), А (К + Ь) = А (К), иначе

— Продолжить вычисления для следующего Ь.

Комментарий. При завершении данного шага в каждом процессе К = 1, 2, 3,..., N значение С (К) будет равно номеру столбцового ведущего процесса

найти столбец

ік ■ %ік = А{К)1ку - с(К)ік < 0;

А (К) = ггк, С (К) = %к

положить С (К) = С (К + Ь), А (К) = А (К + Ь)

если {К : = 0} = 0,

в противном случае,

Ш!аг 3. Если Кс = 0, то сформировать решение задачи:

положить ж[1 : п] = 0;

каждому процессу К = 1, 2, 3,..., N для

положить

хм (К)г = Хв(К)г;

— вернуть х - оптимальное решение задачи, у - оптимальное решение двойственной задачи и завершить выполнение алгоритма.

• Шаг 4. Процессу Кс разослать всем процессам К = 1, 2,..., N столбец Аікс и значение коэффициента С{кс.

• Шаг 5. Каждому процессу К = 1, 2,..., N вычислить д(К) = В^К)-1 Аікс и

г(К) = аге тіп

I: д(К)1>0

Іг(К, I) =

Хв (К), д(К )і

Если г (К) не найдено, то положить С (К) = 0, в противном случае положить С (К) = г (К), А (К) = к(К,г(К)).

• Шаг 6. Для Ь = 1, 2, 4, 8,..., N каждому процессу К = 1,2,3,..., N, номер которого удовлетворяет условию ((К — 1)%(2Ь)) < Ь, осуществить обмен данными с процессом К + Ь:

— Если С (К) = 0, то положить С (К) = С (К + Ь), А (К) = А (К + Ь) и продолжить вычисления для следующего Ь.

— Если С (К + Ь) = 0, то положить С(К + Ь) = С (К), А (К + Ь) = А(К) и нродол-жить вычисления для следующего Ь.

— Если А(К) < А (К + Ь), то положить (7(1^ + Ь) = С (К), А (К + Ь) = А(К), иначе положить С(К) = С(К + Ь), А(К) = А(К + Ь).

— Продолжить вычисления для следующего Ь.

Комментарий. При завершении данного шага в каждом процессе К = 1, 2, 3, значение С(К) будет равно номеру строкового ведущего процесса

Кг = аге1ех тт й,(К,г(К)).

К: Зг(К)

N

• Шаг 7. Ведущему процессу Кт разослать всем процессам К = 1, 2, 3,..., N ведущую строку г* = г(Кт) обратной матрицы

(в{к)-1у ] = (ъг*1, ът,2, •••, ьг*гп )

и значение д(Кг)г*. Положить М(К)Г* = %кс, св{К)г* = С{кс.

• Шаг 8. Каждому процессу К = 1,2,... ,М вычислить новые значения базисных переменных процесса Хв (К) по формулам (5), новый блок обратной матрицы процесса В-1 (К) по формулам (6), и двойственное субрешение блока

у(К) = (В(К)-1)Тсв(К).

• Шаг 9. Для Ь = 1, 2, 4, 8,..., N каждому процессу К = 1,2,3,..., N, номер которого удовлетворяет условию ((К — 1)%(2Ь)) < Ь, осуществить обмен данными с процессом К + Ь:

- у = у(К) + у(К + Ь);

- у(К) = у(К + I) = у-

Таблица 2

Алгебраический вычислительный ресурс метода метода обратной матрицы

Один процессор N процессоров

Количество Количество Нагрузка

Оператор алгебраических пересылаемых на один

операций операндов процесс

Проверка условия оптимальности 2ш[1, п] - 2т[1, п/Х ]

Выбор с-ведущего процесса - N 1о§2 N

Передача ведущего столбца - т + 1 -

Нахождение д ш(ш — 1) - т(т — 1)/Х

Нахождение ведущей строки 2т N 2т/Х + 1о§2 N

Модификация Хв 1 + 2т 1 2т/Х

Модификация В-1 3т2 т 3(т2^)

Вычисление у(К) - - (т/А0(2(т/А0 - 1)

Вычисление у т(2т — 1) N 1о§2 N

Итого: 2т[1, п] + 6 т2 + 2т — 1 2т + 3N + 2 2т[1, п/Х ] + 6т2/N + 2т/Х+ 31о§2 N

— Продолжить вычисления для следующего Ь.

Комментарий. При завершении данного шага в каждом процессе К = 1, 2, 3,..., N вектор у(К) будет содержать двойственное решение у, построенное по текущему базисному решению.

• Конец алгоритма

Таким образом, описание параллельной версии метода обратной матрицы оказывается сложнее описания табличного симплекс-метода. В основном это объясняется большей степенью неоднородности используемых структур. При выполнении одной итерации используется следующие широковещательные коммуникации между процессами: (1) 1о§2 ^] широковещательных коммуникаций между процессами при нахождении столбцового ведущего процесса Кс; (2) широковещательная коммуникация для передачи вводимого в базис (ведущего) столбца А{* и его номера, (3) 1о§2 [-/V] широковещательных коммуникаций между процессами при нахождении строкового ведущего процесса Кг; (4) широковещательная

коммуникация для передачи ведущей строки ^В(К) ^ и ее номера, (5) 1о§2 \М~\ широковещательных коммуникаций между процессами для обмена двойственными субрешениями и нахождения двойственного решения.

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

3. Заключение

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

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

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

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

Работа проводилась при финансовой поддержке РФФИ (проект 10-07-96003-р_урал_а) Статья рекомендована к публикации программным комитетом международной научной конференции «Параллельные вычислительные технологии 2011».

Литература

1. Гаранжа, В.А. Параллельная реализация метода Ньютона для решения больших задач линейного программирования / В.А. Гаранжа, А.И. Голиков, Ю.Г. Евтушенко // Журн.

2. Hall, Ju. Towards a practical parallelization of the simplex method / Ju. Hall // J.

3. Panyukov, A.V. Exact and Guaranteed Accuracy Solutions of Linear Programming Problems by Distributed Computer Systems with MPI / A. V. Panyukov, V. V. Gorbik // Tambov University REPORTS: A Theoretical and Applied Scientific Journal. Series:

http: //vestnik.tsutmb.ru/old/index.php?module=subjects&func=viewpage&pageid=165

4. Панюков, А.В. Библиотека классов «Exact Computational / А.В. Панюков, М.И. Гер-маненко, В.В. Горбик // Программы для ЭВМ, базы данных, топологии интегральных микросхем: официальный бюллетень Рос. агентства по патентам и товарным знакам

5. Панюков, А.В. Параллельные алгоритмы решения систем линейных алгебраических уравнений с применением вычислений без округления / А.В. Панюков, М.И. Герма-ненко, В.В. Горбик // Параллельные вычислительные технологии (ПАВТ—2007): тр.

6. Панюков, А.В. Приложение для безошибочного нахождения обобщенной обратной матрицы методом Мура-Пенроуза и безошибочное решение систем линейных алгебраических уравнений / А.В. Панюков, М.И. Германенко // Информационные технологии мо-

7. Васильев, Ф.П. Линейное программирование / Ф.П. Васильев, А.Ю. Иваницкий - М.:

Анатолий Васильевич Панюков, доктор физико-математических наук, профессор, кафедра «Экономико-математические методы и статистика:», Южно-Уральский государственный университет, [email protected]

Василий Владимирович Горбик, кафедра «Экономико-математические методы и статистика:», Южно-Уральский государственный университет.

Поступила в редакцию 20 марта 2011 г.