§35. Алгоритм решения СЛАУ методом Гаусса

Рассмотрим систему линейных алгебраических уравнений

                                                                                 (1)

с невырожденной матрицей размером . Наиболее известной формой гауссова исключения является та, в которой система (1) приводится к верхнетреугольному виду путем вычитания одних уравнений, умноженных на под­ходящие числа, из других уравнений, полученная треугольная система решается с помощью обратной подстановки. Математи­чески все это эквивалентно тому, что вначале строится разло­жение , где L - нижнетреугольная матрица, с едини­цами на главной диагонали, а U - верхнетреугольная матрица. Затем решаются треугольные системы:

.                                                              (2)

Процесс их решения называется прямой и обратной подста­новками. Рассмотрим сначала LU – разложение, занимающее большую часть времени всего процесса решения СЛАУ. Псевдокод разложения можно представить следующим образом:

Предположим вначале, что мы располагаем вычислительной системой с ло­кальной памятью и числом процессоров p=n. Пусть i-я строка матрицы A хранится в процессоре i. Тогда один из возможных вариантов организации LU-разложения заключается в следующем: на первом шаге первая строка рассылается всем процессорам, после чего вычисления

могут выполняться параллельно процессорами p₂, ...,p_n. На втором шаге вторая строка приведенной матрицы рассылается из процессора p₂ процессорам p₃, ..., p_n , a затем проводятся параллельные вычисления и т. д. Алгоритм первых четырех шагов приведен ниже.

Шаг 1:

Шаг 2: на

Шаг 3:

Шаг 4: на

Отметим два главных недостатка этого подхода: значительный объем обмена данными между процессорами и уменьшение числа активных процессоров на 1 через каждый шаг.

Альтернативой хранению по строкам является вариант, в котором i-й столбец матрицы А хранится в процессоре i. В этом случае на первом шаге все множители l_i1 вычисляются в про­цессоре 1 и рассылаются остальным процессорам. Затем про­цессорами 2,…,n параллельно производятся модификации

.

Вычислив множители 1_i1, процессор 1 прекращает работу; во­обще, с каждым шагом число простаивающих процессоров уве­личивается на единицу. Возникает та же, что и в строчно-ориентированном алгоритме, проблема балансировки нагрузки.

Слоистая схема хранения

В более реалистической ситуации, когда , проблема балансировки нагрузки в известной степени смягчается. Предположим, что и применяется хранение по строкам. По­местим первые k строк матрицы А в память процессора 1, сле­дующие k строк в память процессора 2 и т. д. Этот способ хранения назовем блочной схемой. Снова первая строка рас­сылается из процессора 1 остальным процессорам, а затем выполняются необходимые вычисления, однако теперь это делается блоками по k наборов операций в каждом процессоре. Как и прежде, в ходе приведения все большее число процессоров ста­новятся бездействующими, однако отношение общего времени вычислений к времени обменов и времени простоев является возрастающей функцией от k.

Более привлекательна схема хранения, в которой строки, распределенные в разные процессоры, как бы прослаивают друг друга. Будем по-прежнему считать, что , и пусть строки 1, p+1, 2р+1,... хранятся в процессоре 1, строки 2, p+2, 2p+2, ... - в процессоре 2 и т. д. Такой способ хранения мы будем называть циклической слоистой схемой. Проблема простоя процессоров для этой схемы теряет остроту; например, процессор 1 будет работать почти до самого конца приведения, а именно до тех пор, пока не за­кончится обработка строки (k-1)р+1. Разумеется, некото­рая неравномерность загрузки процессоров сохранится. Так, после первого шага строка 1 станет не нужна, и на следующем шаге процессору 1 придется обрабатывать на одну строчку меньше, чем остальным процессорам. Неравномерность сход­ного типа возникает и в том вполне вероятном случае, когда n не кратно числу процессоров.

Тот же принцип прослаивания можно использовать при хра­нении по столбцам: теперь столбцы 1, р+1, …, (k-1)p+1 закреплены за процессором 1, столбцы 2, р+2,…, (k-1)p+2 - за процессором 2 и т. д. Снова все процессоры (при небольшом дисбалансе нагрузки) будут работать почти до конца приведения.

Опережающая рассылка и опережающее вычисление

Обсудим теперь более подробно LU- разложение с использованием циклической слоистой схемы хранения. На первом шаге аннулируются элементы первого столбца. Если процесс ведется обычным последовательным образом, то вычисляются множители l_i1, модифи­цируются соответствующие строки, а затем начинается второй шаг. В начале этого шага процессор, хранящий вторую строку, должен переслать ее остальным процессорам. Следовательно, возникает задержка на время этой пересылки. Очевидным вы­ходом из положения является немедленная рассылка процессором 2 модифицированной второй строки, как только эта строка приняла окончательный вид. Такую стратегию мы назовем опережающей рассылкой. Если можно совместить вычисления и обмены, то все процессоры будут заняты модификациями первого шага, пока идет рассылка. То же самое будет проис­ходить и на других шагах: на k-м шаге (k+1)-я строка рассылается сразу после того, как окончится ее модификация. Насколько выгодна эта стратегия для конкретной вычислительной системы, зависит от топологии ее межпроцессорных связей.

Решение треугольных систем

После выполне­ния LU-разложения нужно решать треугольные системы уравнений. Мы рассмотрим только верхнетреугольную систему Ux=с. Основными методами снова будут обсуждавшиеся ранее столбцовый алгоритм и алгоритм скалярных произведений. Псевдокод решения верхнетреугольной системы уравнений можно представить следующим образом:

Если считать, что матрица U является результатом ранее выполненного LU-разложения, то способ распределения (хранения) матрицы уже предопределен схемой хранения, использовавшейся при разложении. Например, если использовалась строчная циклическая слоистая схема распределения при разложении, то точно таким же образом построчно будет распределена и матрица U. Предположим, что и правая часть системы распределена послойно по процессорам.  Тогда одна из возможных реализаций строчно-ориентированного алгоритма выглядит так:

Шаг 1: на

Шаг 2:

Шаг 3: на

Шаг 4: на

Шаг 5: и

Шаг 6: на

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

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

В таблице 14 приведены временные затраты на LU-разложение и решение верхнетреугольной системы уравнений. Из данной таблицы видно, что решение верхнетреугольной системы уравнений занимает много меньше 1% от общего времени решения СЛАУ и его алгоритм вполне можно не распараллеливать.

Таблица 14

На рис. 58 приводится листинг параллельной программы для решения СЛАУ столбцово-ориентированным методом Гаусса с выбором ведущего элемента.

Рис. 58. Листинг программы, реализующей параллельный алгоритм решения СЛАУ методом Гаусса