Автором предложен новый метод распараллеливания прогонки на основе использования свойства ассоциативности операций перемножения матриц. В отличие от давно известной версии параллельного алгоритма Стоуна, основанного на приёме сдваивания, новый метод имеет те же характеристики устойчивости, что и последовательная прогонка. Предложена также и блочная модификация метода. 1. Введение В данной статье автор предлагает новые способы решения системы линейных алгебраических уравнений1 с трёхдиагональными матрицами, на той же идейной основе, что и известный метод сдваивания Стоуна, - с использованием ассоциативности матричного умножения. Однако использование последовательных фрагментов позволяет автору применить нормировку, что делает новый метод устойчивым.
Как только у исследователей появились в распоряжении вычислительные устройства с
возможностью параллельной работы, так сразу было предложено несколько параллельных
методов решения трёхдиагональных СЛАУ ([
Все перечисленные алгоритмы имеют логарифмическую (относительно размера задачи) длину критического пути, и придуманы в то время, когда специалисты по распараллеливанию ещё предполагали, что в будущем проблема эффективных пересылок между устройствами будет как-то решена. Эти методы также по-разному загружают доступные вычислительные устройства, что создаёт при их реализации дополнительные проблемы и уменьшает реальную эффективность.
Вышеизложенные проблемы подвигли автора статьи к работе по нахождению альтернативных способов распараллеливания прогонки. Автор преследовал две главных цели: разработать алгоритм, достаточно просто отображаемый на современные архитектуры, и при этом лишённый недостатков описанных выше методов, а также обеспечить его устойчивость. Для этого прежде всего им был изучен самый быстрый из перечисленных методов – метод рекурсивного сдваивания Стоуна, основанный на использовании ассоциативности операции перемножения матриц. 2.2 Основания метода Стоуна и возможные пути его изменения1
Итак, введём следующие обозначения. Пусть задача состоит в решении СЛАУ
где
где
1 В изложении метода Стоуна автор использует собственные обозначения, а не обозначения из статей
самого Стоуна [
Аналогично у Стоуна решается и СЛАУ с матрицей U. При этом выполнение сдваивания на данных этапах решения исходной задачи не вызывает роста вычислительной погрешности. Остаётся рассмотреть, как у Стоуна распараллеливается процесс разложения.
Если применить формулу Бине-Коши (см. например [1,2]) к вычислению ведущего главного минора1 Δk матрицы A порядка k, то получим, что, учитывая (4), (5), (6), 1 У Стоуна – qk, без отметки, что это ведущий главный минор. Доказательство формулы (16) он ведёт непосредственно по формулам вычисления элементов матриц, подобно тому, как это сделано здесь, в формулах (38) – (40), но не в блочном варианте. (22) (23) что схоже с (14) и даёт возможность применить тот же самый приём сдваивания, что и при решении двухдиагональных СЛАУ. Однако у этого этапа есть два существенных отличия. Первое – у произведения нижняя строка остаётся той же, что и у всех сомножителей – , а у произведения нижняя строка та же, что верхняя у . Второе отличие, пожалуй, самое существенное. Дело в том, что если распараллеленные по Стоуну решатели двухдиагональных СЛАУ имеют тот же диапазон устойчивости, что и последовательные, то распараллеленный по Стоуну этап разложения матрицы имеет гораздо худшую устойчивость, чем LU-разложение по компактной схеме метода Гаусса. Дело в том, что, по [1,2], оценки эквивалентного возмущения у нераспараллеленного разложения определяются в основном ростом (или отсутствием роста) элементов полученного разложения. А вот при вычислении элементов матриц рост промежуточных результатов будет гораздо больше. Это и понятно, если вспомнить, что ведущие главные миноры матрицы равны произведениям диагональных элементов получаемых матриц.
Именно поэтому алгоритм Стоуна, несмотря на то, что его часто рассказывают студентам на курсах по параллельным вычислениям, на практике не применяется. Вместе с тем нельзя полностью ставить на нём крест. Нередки случаи, когда после один раз выполненного разложения СЛАУ с уже разложенной матрицей решаются снова и снова, с новой правой частью. Тогда вполне разумным представляется, потратив большое время на это разложение обычным способом, затем использовать элементы метода Стоуна при решении новых потоков СЛАУ.
Кроме того, и эти «куски» метода можно модифицировать так, чтобы они лучше отображались на архитектуру вычислительных комплексов. Как нетрудно видеть по Рис. 1, алгоритм сдваивания имеет граф с довольно сложной структурой передач. Рис. 1. Алгоритм сдваивания для вычисления всех частичных "произведений" для ассоциативных операций при n=8, чёрным обозначены операции, результаты которых нужны на выходе алгоритма Другим важным его недостатком является большой коэффициент избыточности. Для метода сдваивания вычисление всех частных произведений будет стоить вместо n-1 операций умножения матриц таких же операций, так что коэффициент избыточности равен Дополнительную добавку несёт замена операции типа a+bc на, хоть и урезанную благодаря специфике матриц, но всё же занимающую больше ресурсов операцию перемножения матриц. Это накладывает жёсткие требования на количество устройств, нужных для того, чтобы хотя бы не проиграть последовательному алгоритму.
Пришедшая автору в голову идея для замены фрагментов метода Стоуна, связанных с решениями двухдиагональных СЛАУ – замена сдваивания на последовательно-параллельный метод организации вычислений. Как видно на Рис. 12, структура передач при его использовании существенно упрощается. Кроме этого, коэффициент избыточности существенно уменьшается по сравнению со сдваиванием и становится равным 2. Это означает, что данную схему можно применять для ускорения соответствующих частей прогонки уже при небольшом количестве доступных устройств. Рис. 2. Алгоритм последовательно-параллельного метода для вычисления всех частичных "произведений" для ассоциативных операций при n=25, чёрным обозначены операции, результаты которых нужны на выходе алгоритма Есть и ещё важный момент, который следует подчеркнуть. Сам алгоритм сдваивания Стоуна нельзя применять по аналогии, если у нас СЛАУ не с трёхдиагональной, а с блочнотрёхдиагональной матрицей. Однако его фрагменты, связанные с решением двухдиагональных СЛАУ, вполне годятся для ускорения решения СЛАУ с блочно-двухдиагональными матрицами. Применение в последнем случае последовательно-параллельной схемы облегчит адаптацию этих частей к решению многих задач. 3. Новый метод распараллеливания главной части прогонки - разложения матрицы.
Алгоритм сдваивания Стоуна неустойчив из-за возможного роста результатов при выполнении промежуточных вычислений, связанных с ведущими главными минорами матрицы. Однако нам не нужны сами эти миноры, нужны лишь отношения соседних. Посмотрим снова на формулы, связывающие их, – нельзя ли как-то избежать роста результатов промежуточных вычислений? 3.1 Идея введения нормировки в промежуточных вычислениях
Посмотрим на формулы (21) – (23) снова. Выполним в (21) подстановку не до конца, как в (23), а до некоторого значения i<k: Обозначим произведения : (24) Обозначим элементы верхней строки этой матрицы как видно, что ( , (26) (27) (27')
1 (29') 2 и таким образом мы уже избавились от одного из источников больших промежуточных результатов – самих миноров. Остаётся избавиться от вычислений с большими элементами промежуточных матриц и . Если посмотреть на (28) внимательно, то окажется, что вместо 2мерных строк и в ней могут фигурировать эти же строки, домноженные на любой выбранный нами нормировочный коэффициент. Вкупе с идеей последовательно-параллельного метода это даёт следующую схему. 3.2 Схема нового метода разложения трёхдиагональной матрицы
Весь диапазон натуральных чисел от 1 до n разбивается на q промежутков – от 1 до k1, от k1+1 до k2, ..., от kq-1+1 до kq=n. На первом промежутке все значения ui вычисляются последовательно, по стандартным формулам разложения на 2 диагональные матрицы. На остальных промежутках выполняется следующая процедура.
Пусть рассматривается j+1-й промежуток. Тогда уже известны значения , , (29)
, Формулы (30) – (31') как раз включают в себя нормировку – их выполнение приводит к единице один из коэффициентов в знаменателе дроби в (28). После того, как выполнены вычисления на всех промежутках, на них (кроме первого) ещё неизвестны значения ui. Их мы для каждого j+1-го промежутка вычисляем параллельно по формулам
Естественно, что каждое значение можно вычислить только после вычисления . В результате, если разбиение на интервалы будет проведено равномерно, структура алгоритма разложения будет иметь почти тот же вид, что и на Рис. 12, с тем, однако, отличием, что «тяжести» операций будут разными. Это отражено на Рис. 13.
Разнородность операций делает не столь тривиальной задачу разбиения выполняемых операций по разным устройствам; тем не менее, эта задача существенно проще, чем организация пересылок в методах, использующих сдваивание. Оценим количество операций на каждом этапе. Вычисления будем считать выполненными заранее. 1 последний переход получается делением числителя и знаменателя на 2 деление на себя не показано 1. Вычисление разложения на первом из промежутков известно, если длина промежутка равна m, то это по m-1 операций деления, умножения и сложения/вычитания. Все эти операции следуют друг за другом последовательно.
2. Вычисления коэффициентов на остальных промежутках, если их длина равна m, займёт m-2 параллельных операций: одна из них – умножение, деление, потом вычитание/сложение, вторая – два параллельных умножения, вычитание/сложение, потом деление, третья – деление.
Видно, что вычисление коэффициентов в общем сложнее. В него можно добавить и вычисление числителя последней дроби из (32), он не зависит от значений . Тогда для вычислений на 2м этапе самым длинным является вычисление по «низу картинки» – это операция типа сложение, потом деление и потом снова сложение, повторяемая q-1 раз. Параллельно этой операции нужно выполнить ещё m-2 таких же, а также ещё одну, состоящую только из деления и сложения.
Приведённые формулы опираются на формулы (19) для миноров и поэтому, как и сам метод сдваивания Стоуна, казалось бы, неприменимы для блочно-трёхдиагональных матриц. Рис. 3. Алгоритм нового метода для вычисления разложения трёхдиагональной матрицы при n=25 и равномерном разбиении на 5 промежутков, разная интенсивность соответствует разным операциям 4. Блочная версия
Пусть теперь нам нужно выполнить разложение блочно-трёхдиагональной матрицы где где , , – квадратные блоки одинакового размера. Для разложения (34) существуют формулы, являющиеся частью блочной прогонки: а после домножения этого равенства справа на произведение это даст нам , (40) Теперь наличие двучленной рекурсии позволяет нам повторить приём объединения блоков в блочные «вектора»: и получить аналогичную (12) формулу , (35) (36) (38) (39) (41) (42) (43) (44) (45) вид(46) (47) (47') где где Теперь мы можем повторить рассуждения, аналогичные формулам из 3.1. Выполним в (42) подстановки не до конца, а до некоторого значения i<k: обозначим произведения
: обозначим блоки верхней «строки» этой матрицы как но, что и . Тогда из вида матрицы ( ,
). Поэтому Учитывая (40), получаем (48) (49) (50) (51) (52) (53) (54) и, таким образом, нам не нужно вычислять произведения матриц. От роста же матриц и можно попытаться избавиться с помощью нормировки. К сожалению, она не может быть вполне подобна нормировке в (3.2), и вопрос об её устойчивости остаётся открытым. Приведём возможную схему.
Весь диапазон натуральных чисел от 1 до n разбивается на q промежутков – от 1 до k1, от k1+1 до k2, ..., от kq-1+1 до kq=n. На первом промежутке все значения Ui вычисляются последовательно, по стандартным формулам (37). На остальных промежутках выполняется следующая процедура.
Пусть рассматривается j+1-й промежуток. Тогда уже известны значения , , (55) , (55’) удовлетворяет заявленному требованию скалярности наддиагональных блоков и потому может быть разложена предложенным методом. К сожалению, ограничение невырожденности для блоков, составляющих хотя бы одну из побочных диагоналей, всё же довольно сильно ограничивает область применимости блочной версии последовательно-параллельного метода разложения блочно-трёхдиагональной матрицы.
После этого вместо систем вида (1) можно будет либо решать СЛАУ вида Формулы (30) – (31') как раз включают в себя нормировку. После того, как выполнены вычисления на всех промежутках, на них (кроме первого) ещё неизвестны значения Ui. Их мы для каждого j+1-го промежутка вычисляем параллельно по формулам
Естественно, что каждое значение можно вычислить только после вычисления . В результате, если разбиение на интервалы будет проведено равномерно, структура алгоритма разложения будет иметь тот же вид, что и на Рис. 13. Естественно, на этот раз вершинам графа будут соответствовать более сложные операции.
Конечно, скалярность блоков – довольно сильное ограничение для того, чтобы можно было распараллелить блочную прогонку. На деле можно ограничиться менее сильным ограничением – их невырожденностью. Действительно, пусть все наддиагональные блоки матрицы невырождены. Тогда вместо матрицы A мы можем выполнить предложенным способом разложение матрицы DA, где Действительно, матрица (59) (60) (61) (62) (63) либо домножить матрицу L из полученного LU-разложения матрицы DA на матрицу D-1 слева и работать с полученным LU-разложением.
Остаётся отметить, что решения СЛАУ с полученными из разложения блочнодвухдиагональными матрицами распараллеливаются по схеме, принципиально не отличающейся от схемы, получающейся из формул (10) – (14), с заменой скалярных операций на блочные, и с возможностью упрощения по последовательно-параллельной схеме, как на Рис. 12. В отличие от блочной схемы разложения, при этом распараллеливании не нужно применять специальные методы для обеспечения устойчивости – она у распараллеленной схемы та же, что и у последовательных блочных версий. Автор при разработке последовательно-параллельного метода для разложения трёхдиагональной матрицы, изложенного в части 3, прежде всего руководствовался целью распараллелить нахождение того же LU-разложения трёхдиагональной матрицы, что получается последовательной схемой, получаемой из компактной схемы метода Гаусса применительно к трёхдиагональным матрицам. В условиях точных вычислений предложенный метод эквивалентен как указанной версии компактной схемы метода Гаусса, так и методу рекурсивного сдваивания Стоуна. Структура расположения ненулевых элементов матрицы A дублируется аналогичной структурой разложения (см. Рис. 14). Рис. 4. Структура расположения ненулевых элементов исходной матрицы. В нижнем треугольнике она дублируется структурой матрицы L, в верхнем – структурой матрицы U. Другие алгоритмы, возможно, похожие по схеме, дают другие разложения. Например, один из читателей увидел схожесть предложенного метода и следующего: «перестановки такие, что в трехдиагональной матрице разделители уходят в конец матрицы; возникает набор не связанных блоков и окаймление; такую матрицу можно устойчиво факторизовать в параллельном режиме» (возможный пример такого с одним разделением приведён на Рис. 15). Рис. 5. Выполнение перестановок в окаймлении, ассоциации с которым могут возникнуть у читателя статьи. Выполнено только одно разделение матрицы на фрагменты На Рис. 16 показана структура ненулевых элементов получаемого разложения. Рис. 6. Структура расположения ненулевых элементов разложения «переставленной» матрицы. В нижнем треугольнике она дублируется структурой матрицы L, в верхнем – структурой матрицы U. Хорошо видно появление дополнительных ненулевых элементов. При разбиении перестановками верхнего левого трёхдиагонального блока исходной матрицы будут возникать аналогичные структуры. Кроме этого, таким методом, в отличие от авторского, вычисляется LUразложение не для исходной матрицы A, а для матрицы PAPT, где P – матрица перестановок .Поэтому такая схема явно неэквивалентна предложенному автором методу, в котором вычисляемое разложение не имеет таких дополнительных элементов. При этом, однако, подмеченная читателем схожесть в том, что основными частями разложения и там, и там является параллельная обработка отдельных трёхдиагональных фрагментов исходной матрицы, плюс некоторое добавление. Однако и обработка фрагментов, и добавления – разные в разных методах. К слову, несмотря на дополнительное «заплывание» структуры разложения в методе такого окаймления тот, пожалуй, более пригоден для распараллеливания разложения блочнотрёхдиагональной матрицы в тех случаях, когда у нас нет невырожденности блоков хотя бы на одной из её побочных диагоналей. 6. Заключение
Переход от приёма сдваивания к варианту последовательно-параллельного выполнения для использования ассоциативности операций позволил сконструировать такой метод, где, благодаря использованию приёма нормировки, характерного для последовательных операций [2], возможно стало сочетать параллельность с устойчивостью разложения трёхдиагональной матрицы. Аналогичный приём, который автор предлагает использовать и при решении двухдиагональных СЛАУ, получающихся после разложения, также позволяет упростить схему организации вычислений и снять жёсткие требования на количество устройств, позволяя запрограммировать эффективное параллельное исполнение даже при сравнительно небольшом уровне параллелизма вычислительной системы.
У читателей, возможно, возникнут вопросы по поводу апробации метода и его сравнению с другими методами решения этой же задачи. На апробацию у автора просто пока не было времени (контуры метода намечены в апреле, а в нынешнем виде он записан в середине мая 2015г.), но в ближайшем будущем запланировано заняться и апробацией метода, и детальным его сравнением с другими.
Кроме получения непосредственно новых схем вычисления, автор рекомендует читателю присмотреться к последовательно-параллельной схеме вычислений там, где схемы сдваивания не дают возможности построить устойчивые алгоритмы. Не исключено, что подобные замены могут помочь и в других аналогичных случаях. Литература Yet another tridiagonal matrix algorithm parallelizing method Alexey Frolov Keywords: Thomas algorithm, tridiagonal matrix algorithm, parallelizing New tridiagonal matrix algorithm parallelizing method is described in this article.