Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org Ещё один метод распараллеливания прогонки с использованием ассоциативности операций* А.В. Фролов ИВМ РАН Автором предложен новый метод распараллеливания прогонки на основе использо- вания свойства ассоциативности операций перемножения матриц. В отличие от давно известной версии параллельного алгоритма Стоуна, основанного на приёме сдваива- ния, новый метод имеет те же характеристики устойчивости, что и последовательная прогонка. Предложена также и блочная модификация метода. 1. Введение В данной статье автор предлагает новые способы решения системы линейных алгебраиче- ских уравнений1 с трёхдиагональными матрицами, на той же идейной основе, что и известный метод сдваивания Стоуна, – с использованием ассоциативности матричного умножения. Одна- ко использование последовательных фрагментов позволяет автору применить нормировку, что делает новый метод устойчивым. 2. Прогонка как метод решения СЛАУ специального вида. Прогонка – последовательный алгоритм решения трёхдиагональной СЛАУ – является ча- стным случаем общего метода исключения неизвестных, однако получила специальное назва- ние из-за распространённости задач такого типа в прикладных исследованиях. В виде отдель- ного алгоритма «открыта» несколькими исследователями независимо друг от друга (И.М. Гельфандом и О.В. Локуциевским в СССР, Л.Х. Томасом в США2). С этими названиями описа- на в большом количестве учебников по численным методам (например, в [5]), причём в разных вариантах. Классическая схема прогонки практически не содержит возможностей для распараллели- вания, максимум – по двум потокам вычислений. Однако решение трёхдиагональных СЛАУ довольно востребовано в различных вычислительных задачах при моделировании, поэтому за- дача «распараллеливания прогонки», т.е. параллельного решения исходной задачи, давно при- влекает внимание многих исследователей. В подавляющем большинстве учебников (например, в [5]) отмечается, что задача решения СЛАУ Ax=b с трёхдиагональной матрицей A при решении с помощью прогонки эквивалентна последовательности двух задач. Это разложение матрицы в произведение, например, двухдиа- гональных матриц, а также решение СЛАУ с этими матрицами. Такая разбивка впоследствии даёт возможность использовать уже найденное разложение для более быстрого решения СЛАУ и с другими правыми частями. В дальнейшем мы увидим это как на примере рассмотрения ал- горитма Стоуна, так и при конструировании собственного алгоритма. 2.1 Известные параллельные методы решения трёхдиагональных СЛАУ Как только у исследователей появились в распоряжении вычислительные устройства с возможностью параллельной работы, так сразу было предложено несколько параллельных ме- тодов решения трёхдиагональных СЛАУ ([7–9]). Сравнение с позиций середины 70-х гг. XX * Исследование выполнено при частичной финансовой поддержке гранта Российского научного фонда (проект N14-11-00190). 1 далее для систем линейных алгебраических уравнений будем использовать сокращение СЛАУ. 2 поэтому в англоязычной литературе, кроме названия «Tridiagonal matrix algorithm», используется и на- звание «Thomas algorithm» 151 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org века трёх таких методов – циклической редукции, Бунемана и рекурсивного сдваивания – мож- но прочитать в работе автора последнего Х.Стоуна 1975г. [10]. Несмотря на общую привязку этого сравнения к архитектурам того времени, эти сравнения актуальны и теперь. Согласно им, алгоритм рекурсивного сдваивания (Стоуна) имеет лучшие характеристики по количеству тре- буемых вычислительных затрат. Однако читателям, должно быть, известно, что на практике применяют не его, а метод циклической редукции. Это связано с тем, что один из этапов метода рекурсивного сдваивания Стоуна (разложение матрицы в произведение двухдиагональных мат- риц) имеет гораздо худшую устойчивость. Все перечисленные алгоритмы имеют логарифмическую (относительно размера задачи) длину критического пути, и придуманы в то время, когда специалисты по распараллеливанию ещё предполагали, что в будущем проблема эффективных пересылок между устройствами бу- дет как-то решена. Эти методы также по-разному загружают доступные вычислительные уст- ройства, что создаёт при их реализации дополнительные проблемы и уменьшает реальную эф- фективность. Вышеизложенные проблемы подвигли автора статьи к работе по нахождению альтернатив- ных способов распараллеливания прогонки. Автор преследовал две главных цели: разработать алгоритм, достаточно просто отображаемый на современные архитектуры, и при этом лишён- ный недостатков описанных выше методов, а также обеспечить его устойчивость. Для этого прежде всего им был изучен самый быстрый из перечисленных методов – метод рекурсивного сдваивания Стоуна, основанный на использовании ассоциативности операции перемножения матриц. 2.2 Основания метода Стоуна и возможные пути его изменения1 Итак, введём следующие обозначения. Пусть задача состоит в решении СЛАУ (1) где (2) – трёхдиагональная матрица, для которой выполняются условия устойчивости решения мето- дом прогонки [5], (3) – вектор правой части. Если теперь матрицу A разложить в произведение двух треугольных (4) где (5) 1 В изложении метода Стоуна автор использует собственные обозначения, а не обозначения из статей самого Стоуна [9, 10] 152 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org (6) то (1) будет заменена на (7) После этого сначала нужно решить систему (8) и после неё (9) Занявшись сначала более простым распараллеливанием решения двухдиагональных СЛАУ, видим, что для (8) получаются формулы , , (10) после чего введением векторов (11) достигается выражение , (12) где (13) После подстановки в (12) этой же формулы с меньшими значениями индекса получаем (14) после чего нужно вычислить все такие частные произведения, что делается, например, методом сдваивания за ярусов. Аналогично у Стоуна решается и СЛАУ с матрицей U. При этом выполнение сдваивания на данных этапах решения исходной задачи не вызывает роста вычислительной погрешности. Остаётся рассмотреть, как у Стоуна распараллеливается процесс разложения. Если применить формулу Бине-Коши (см. например [1,2]) к вычислению ведущего главно- го минора1 Δk матрицы A порядка k, то получим, что, учитывая (4), (5), (6), и, вводя Δ0 = 1, имеем (16) а из формул перемножения матриц - (17) и (18) Остаётся понять, как параллельно вычислить все ведущие главные миноры матрицы A. Из учебников и справочников (например, [1,2]) известно, что у трёхдиагональных матриц для ве- дущих главных миноров выполняются рекуррентные соотношения (19) что позволяет, введя вектора (20) получить аналогичную (12) формулу , (21) где 1 У Стоуна – qk, без отметки, что это ведущий главный минор. Доказательство формулы (16) он ведёт непосредственно по формулам вычисления элементов матриц, подобно тому, как это сделано здесь, в формулах (38) – (40), но не в блочном варианте. 153 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org (22) и (23) что схоже с (14) и даёт возможность применить тот же самый приём сдваивания, что и при ре- шении двухдиагональных СЛАУ. Однако у этого этапа есть два существенных отличия. Пер- вое – у произведения нижняя строка остаётся той же, что и у всех сомножителей – , а у произведения нижняя строка та же, что верхняя у . Второе отличие, пожалуй, самое существенное. Дело в том, что если распараллеленные по Стоуну ре- шатели двухдиагональных СЛАУ имеют тот же диапазон устойчивости, что и последователь- ные, то распараллеленный по Стоуну этап разложения матрицы имеет гораздо худшую устой- чивость, чем LU-разложение по компактной схеме метода Гаусса. Дело в том, что, по [1,2], оценки эквивалентного возмущения у нераспараллеленного разложения определяются в основ- ном ростом (или отсутствием роста) элементов полученного разложения. А вот при вычисле- нии элементов матриц рост промежуточных результатов будет гораздо больше. Это и понятно, если вспомнить, что ведущие главные миноры матрицы равны произведениям диа- гональных элементов получаемых матриц. Именно поэтому алгоритм Стоуна, несмотря на то, что его часто рассказывают студентам на курсах по параллельным вычислениям, на практике не применяется. Вместе с тем нельзя полностью ставить на нём крест. Нередки случаи, когда после один раз выполненного разложе- ния СЛАУ с уже разложенной матрицей решаются снова и снова, с новой правой частью. Тогда вполне разумным представляется, потратив большое время на это разложение обычным спосо- бом, затем использовать элементы метода Стоуна при решении новых потоков СЛАУ. Кроме того, и эти «куски» метода можно модифицировать так, чтобы они лучше отобража- лись на архитектуру вычислительных комплексов. Как нетрудно видеть по Рис. 1, алгоритм сдваивания имеет граф с довольно сложной структурой передач. Рис. 1. Алгоритм сдваивания для вычисления всех частичных "произведений" для ассоциативных опе- раций при n=8, чёрным обозначены операции, результаты которых нужны на выходе алгоритма Другим важным его недостатком является большой коэффициент избыточности. Для ме- тода сдваивания вычисление всех частных произведений будет стоить вместо n-1 операций ум- ножения матриц таких же операций, так что коэффициент избыточности равен . Дополнительную добавку несёт замена операции типа a+bc на, хоть и урезанную благодаря 154 Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org специфике матриц, но всё же занимающую больше ресурсов операцию перемножения матриц. Это накладывает жёсткие требования на количество устройств, нужных для того, чтобы хотя бы не проиграть последовательному алгоритму. Пришедшая автору в голову идея для замены фрагментов метода Стоуна, связанных с ре- шениями двухдиагональных СЛАУ – замена сдваивания на последовательно-параллельный ме- тод организации вычислений. Как видно на Рис. 12, структура передач при его использовании существенно упрощается. Кроме этого, коэффициент избыточности существенно уменьшает- ся по сравнению со сдваиванием и становится равным 2. Это означает, что данную схему мож- но применять для ускорения соответствующих частей прогонки уже при небольшом количестве доступных устройств. Рис. 2. Алгоритм последовательно-параллельного метода для вычисления всех частичных "произведе- ний" для ассоциативных операций при n=25, чёрным обозначены операции, результаты которых нужны на выходе алгоритма Есть и ещё важный момент, который следует подчеркнуть. Сам алгоритм сдваивания Сто- уна нельзя применять по аналогии, если у нас СЛАУ не с трёхдиагональной, а с блочно- трёхдиагональной матрицей. Однако его фрагменты, связанные с решением двухдиагональных СЛАУ, вполне годятся для ускорения решения СЛАУ с блочно-двухдиагональными матрица- ми. Применение в последнем случае последовательно-параллельной схемы облегчит адаптацию этих частей к решению многих задач. 3. Новый метод распараллеливания главной части прогонки - разло- жения матрицы. Алгоритм сдваивания Стоуна неустойчив из-за возможного роста результатов при выпол- нении промежуточных вычислений, связанных с ведущими главными минорами матрицы. Од- нако нам не нужны сами эти миноры, нужны лишь отношения соседних. Посмотрим снова на формулы, связывающие их, – нельзя ли как-то избежать роста результатов промежуточных вы- числений? 3.1 Идея введения нормировки в промежуточных вычислениях Посмотрим на формулы (21) – (23) снова. Выполним в (21) подстановку не до конца, как в (23), а до некоторого значения i