Обратная подстановка (вещественный вариант): различия между версиями
| [непроверенная версия] | [досмотренная версия] |
Frolov (обсуждение | вклад) |
ASA (обсуждение | вклад) Нет описания правки |
||
| (не показано 60 промежуточных версий 7 участников) | |||
| Строка 1: | Строка 1: | ||
== | {{algorithm | ||
| name = Обратная подстановка для неособенной верхней треугольной матрицы | |||
| serial_complexity = <math>O(n^2)</math> | |||
| pf_height = <math>O(n)</math> | |||
| pf_width = <math>O(n)</math> | |||
| input_data = <math>O(n^2)</math> | |||
| output_data = <math>n</math> | |||
}} | |||
Основные авторы описания: [[Участник:Frolov|А.В.Фролов]]. | |||
== Свойства и структура алгоритма == | |||
=== Общее описание алгоритма === | |||
= | '''Обратная подстановка''' - решение ''системы линейных алгебраических уравнений'' ('''СЛАУ''') <math>Ux = y</math> с верхней треугольной матрицей <math>U</math>. Матрица <math>U</math> может быть одной из составляющих матрицы <math>A</math> в каких-либо разложениях и получается либо из <math>LU</math>-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, [[Метод Холецкого (квадратного корня), точечный вещественный вариант|разложение Холецкого]] и др.), либо из других (например из QR-разложения). В силу треугольности <math>U</math> решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами. | ||
Исходные данные: | В<ref>В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. - М.: Наука, 1984.</ref> методом '''обратной подстановки''' назван также и метод решения СЛАУ с ''нижней треугольной матрицей''. Там же отмечено, что в литературе иногда под ''обратной подстановкой'' имеют в виду, как и здесь, только решения СЛАУ с ''верхней треугольной матрицей'', а решение ''нижних'' треугольных систем называют [[Прямая подстановка (вещественный вариант)|прямой подстановкой]]. Такой же системы названий будем придерживаться и здесь, во избежание одноимённого названия различных алгоритмов. Кроме того, '''обратная подстановка''', представленная здесь, одновременно может быть частью '''метода Гаусса для решения СЛАУ''', а именно - его '''обратным ходом''', чего нельзя сказать про [[Прямая подстановка (вещественный вариант)|прямую подстановку]]. | ||
Существует метод со сходным названием - [[Обратная подстановка с нормировкой]]. При том, что он решает, по существу, ту же задачу, что и простая '''обратная подстановка''', его схема несколько сложнее. Это связано со специальными мерами по уменьшению влияния ошибок округления на результат. Алгоритм обратной подстановки с нормировкой в данном разделе не рассматривается. | |||
=== Математическое описание алгоритма === | |||
Исходные данные: верхняя треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>). | |||
Вычисляемые данные: вектор решения <math>x</math> (элементы <math>x_{i}</math>). | Вычисляемые данные: вектор решения <math>x</math> (элементы <math>x_{i}</math>). | ||
| Строка 25: | Строка 38: | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
Вычислительное ядро обратной подстановки можно составить из множественных (всего их <math>n-1</math>) вычислений скалярных произведений строк матрицы <math>U</math> на уже вычисленную часть вектора <math>x</math>: | Вычислительное ядро алгоритма обратной подстановки можно составить из множественных (всего их <math>n-1</math>) вычислений скалярных произведений строк матрицы <math>U</math> (за исключением диагонального элемента) на уже вычисленную часть вектора <math>x</math>: | ||
:<math> \sum_{j = i+1}^{n} u_{ij} x_{j} </math> | :<math> \sum_{j = i+1}^{n} u_{ij} x_{j} </math> | ||
в режиме накопления или | в режиме накопления или, в зависимости от требований задачи, без его использования, с их последующим вычитанием из компоненты вектора <math>y</math> и деления на диагональный элемент матрицы <math>U</math>. В отечественных реализациях, даже в последовательных версиях, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм вида | ||
:<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math> | :<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math> | ||
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из | в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента вектора», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому вычислительным ядром алгоритма следует считать не вычисления скалярных произведений, а вычисления выражений | ||
:<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math> | :<math> y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math> | ||
в режиме накопления или | в режиме накопления или, в зависимости от требований задачи, без его использования, плюс деления результатов этих вычислений на диагональные элементы матрицы. | ||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
Как уже | Как уже упоминалось в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть метода обратной подстановки составляют множественные (всего их <math>n-1</math>) вычисления сумм | ||
:<math>y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math> | :<math>y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math> | ||
в режиме накопления или без | в режиме накопления или без его использования, а также деления результатов этих вычислений на диагональные элементы матрицы. | ||
=== | === Схема реализации последовательного алгоритма === | ||
Последовательность действий описываемого алгоритма: | |||
1. <math>x_{n} = y_{n}/u_{nn}</math> | 1. <math>x_{n} = y_{n}/u_{nn}</math> | ||
Далее для всех <math>i</math> от <math>n-1</math> до <math>1</math> по убыванию выполняются | Далее для всех <math>i</math> от <math>n-1</math> до <math>1</math> ''по убыванию'' выполняются | ||
2. <math>x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}</math> | 2. <math>x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}</math>. | ||
Особо отметим, что вычисления сумм вида <math>y_{i} - \ | Особо отметим, что вычисления сумм вида <math>y_{i} - \sum{}_{j = i+1}^{n} u_{ij} x_{j}</math> производят в режиме накопления вычитанием из <math>y_{i}</math> произведений <math>u_{ij} x_{j}</math> для <math>j</math> от <math>n</math> до <math>i + 1</math>, '''''c убыванием''''' <math>j</math>. '''''Использование другого порядка выполнения суммирования приводит к резкому ухудшению параллельных свойств алгоритма''''', хотя, к сожалению, это кое-где встречается в литературе и пакетах программ. В качестве примера использования другого порядка вычислений можно привести фрагмент программы из<ref>Дж.Форсайт, К.Моллер. Численное решение систем линейных алгебраических уравнений. - М.: Мир, 1969.</ref>, где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями используемого авторами книги старого диалекта Фортрана. | ||
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === | ||
Для обратной подстановки порядка n | Для алгоритма обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка <math>n</math> в последовательном (наиболее быстром) варианте требуется: | ||
* <math>n</math> делений, | * <math>n</math> делений, | ||
* | * <math>\frac{n^2-n}{2}</math> сложений (вычитаний), | ||
* <math>\frac{n^2-n}{2}</math> умножений. | |||
Выполнение умножений и сложений (вычитаний) составляет ''основную часть алгоритма''. | |||
При этом использование режима накопления требует | При этом использование режима накопления требует выполнения умножений и вычитаний в режиме двойной точности (или использования функции аналогичной DPROD на языке Фортран). | ||
Таким образом, при классификации по последовательной сложности, метод обратной подстановки относится к алгоритмам ''со сложностью'' <math>O(n^2)</math>. | |||
=== Информационный граф === | === Информационный граф === | ||
| Строка 76: | Строка 91: | ||
Граф алгоритма обратной подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности. | Граф алгоритма обратной подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности. | ||
[[Файл:DirectU.png|450px|thumb|left|Обратная подстановка]] | [[Файл:DirectU.png|450px|thumb|left|Рис. 1. Обратная подстановка]] | ||
'''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция | '''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция выполняет функцию деления. | ||
Естественно введённая единственная координата каждой из вершин <math>i</math> | Естественно введённая единственная координата каждой из вершин <math>i</math> изменяется в диапазоне от <math>n</math> до <math>1</math>, принимая все целочисленные значения. | ||
Делимое в этой операции: | Делимое в этой операции: | ||
| Строка 88: | Строка 103: | ||
Делитель для этой операции - элемент ''входных данных'', а именно <math>u_{nn}</math>. | Делитель для этой операции - элемент ''входных данных'', а именно <math>u_{nn}</math>. | ||
Результат | Результат выполнения операции является ''выходным данным'' <math>x_{i}</math>. | ||
'''Вторая''' группа вершин расположена в двумерной области, соответствующая ей операция | '''Вторая''' группа вершин расположена в двумерной области, соответствующая ей операция <math>a-bc</math>. | ||
Естественно введённые координаты области | Естественно введённые координаты области: | ||
* <math>i</math> — меняется в диапазоне от <math>n-1</math> до <math>1</math>, принимая все целочисленные значения; | * <math>i</math> — меняется в диапазоне от <math>n-1</math> до <math>1</math>, принимая все целочисленные значения; | ||
* <math>j</math> — меняется в диапазоне от <math>n</math> до <math>i+1</math>, принимая все целочисленные значения. | * <math>j</math> — меняется в диапазоне от <math>n</math> до <math>i+1</math>, принимая все целочисленные значения. | ||
Аргументы операции | Аргументы операции: | ||
*<math>a</math>: | *<math>a</math>: | ||
** при <math>j = n</math> элемент входных данных <math>y_{i}</math>; | ** при <math>j = n</math> элемент ''входных данных'' <math>y_{i}</math>; | ||
** при <math>j < n</math> — результат | ** при <math>j < n</math> — результат выполнения операции, соответствующей вершине из второй группы, с координатами <math>i, j+1</math>; | ||
*<math>b</math> — элемент ''входных данных'', а именно <math>u_{ij}</math>; | *<math>b</math> — элемент ''входных данных'', а именно <math>u_{ij}</math>; | ||
*<math>c</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатой <math>j</math> | *<math>c</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатой <math>j</math>. | ||
Результат срабатывания операции является ''промежуточным данным'' алгоритма. | Результат срабатывания операции является ''промежуточным данным'' алгоритма. | ||
Описанный граф | Описанный граф представлен на рис. 1, выполненном для случая <math>n = 5</math>. Здесь вершины первой группы обозначены жёлтым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена только подача входных данных из вектора <math>y</math>, подача элементов матрицы <math>U</math>, идущая во все вершины, на рисунке не представлена. | ||
=== | === Ресурс параллелизма алгоритма === | ||
Для обратной подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы: | Для обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка <math>n</math> в параллельном варианте требуется последовательно выполнить следующие ярусы: | ||
* <math>n</math> ярусов делений (в каждом из ярусов одно деление), | * <math>n</math> ярусов делений (в каждом из ярусов одно деление), | ||
* по <math>n - 1</math> ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от <math>1</math> до <math>n-1</math>. | * по <math>n - 1</math> ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от <math>1</math> до <math>n-1</math>. | ||
| Строка 116: | Строка 131: | ||
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти. | При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти. | ||
Таким образом, при классификации по высоте ЯПФ, алгоритм обратной подстановки относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''линейной''. | |||
=== | === Входные и выходные данные алгоритма === | ||
'''Входные данные''': верхняя треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>). | '''Входные данные''': верхняя треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>). | ||
'''Объём входных данных''': <math>\frac{n (n + 3)}{2}</math> (в силу треугольности достаточно хранить только | '''Объём входных данных''': <math>\frac{n (n + 3)}{2}</math> (в силу треугольности достаточно хранить только ненулевые элементы матрицы <math>U</math>). | ||
'''Выходные данные''': вектор решения <math>x</math> (элементы <math>x_{i}</math>). | '''Выходные данные''': вектор решения <math>x</math> (элементы <math>x_{i}</math>). | ||
| Строка 134: | Строка 149: | ||
При этом вычислительная мощность алгоритма обратной подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь ''константа''. | При этом вычислительная мощность алгоритма обратной подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь ''константа''. | ||
Представленный алгоритм обратной подстановки является полностью детерминированным. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и меняет сложность с линейной на квадратичную. | |||
Наличие линейного количества ярусов ЯПФ, состоящих из одного-единственного деления, потенциально замедляющее параллельные реализации алгоритма, является его характерным "узким местом", особенно в сравнении со схожей по решаемой математической задаче [[Прямая подстановка (вещественный вариант)|прямой подстановке]], где диагональные элементы единичны. В связи с этим для решения СЛАУ предпочтительны такие разложения, содержащие треугольные матрицы, где в треугольных матрицах диагональные элементы единичны. В тех же случаях, когда получаются неособенные треугольные матрицы, их желательно предварительно, до решения СЛАУ с ними, преобразовать в произведение диагональной и треугольной с единичными диагональными элементами. | Наличие линейного количества ярусов ЯПФ, состоящих из одного-единственного деления, потенциально замедляющее параллельные реализации алгоритма, является его характерным "узким местом", особенно в сравнении со схожей по решаемой математической задаче [[Прямая подстановка (вещественный вариант)|прямой подстановке]], где диагональные элементы единичны. В связи с этим для решения СЛАУ предпочтительны такие разложения, содержащие треугольные матрицы, где в треугольных матрицах диагональные элементы единичны. В тех же случаях, когда получаются неособенные треугольные матрицы, их желательно предварительно, до решения СЛАУ с ними, преобразовать в произведение диагональной и треугольной с единичными диагональными элементами. | ||
== Программная реализация == | У алгоритма обратной подстановки существует несколько блочных вариантов. Граф некоторых из них совпадает с графом точечного варианта, различия связаны в основном с порядком прохождения основных циклов алгоритма, а именно, с их развёртыванием и перестановкой. Эти приёмы могут помочь в оптимизации обменов на конкретных вычислительных системах. | ||
== Программная реализация алгоритма == | |||
=== Особенности реализации последовательного алгоритма === | === Особенности реализации последовательного алгоритма === | ||
| Строка 145: | Строка 162: | ||
<source lang="fortran"> | <source lang="fortran"> | ||
X(N) = Y(N) / U (N, N) | X(N) = Y(N) / U(N,N) | ||
DO I = N-1, 1, -1 | |||
S = Y(I) | |||
DO J = N, I+1, -1 | |||
S = S - DPROD(U(I,J), X(J)) | |||
END DO | |||
X(I) = S / U(I,I) | |||
END DO | END DO | ||
</source> | </source> | ||
При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности. | При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности. | ||
=== | === Возможные способы и особенности параллельной реализации алгоритма === | ||
Вариантов параллельной реализации алгоритма не так уж и много, если не использовать то, что оба главных цикла можно развернуть, перейдя, таким образом, к блочной версии. Версии без развёртывания циклов возможны как с полностью параллельными циклами по I: | |||
<source lang="fortran"> | |||
DO PARALLEL I = 1, N | |||
X(I) = Y(I) | |||
END DO | |||
DO J = N, 1, -1 | |||
X(J) = X(J) / U(J,J) | |||
DO PARALLEL I = 1, J-1 | |||
X(I) = X(I) - U(I,J)*X(J) | |||
END DO | |||
END DO | |||
</source> | |||
так и с использованием "скошенного параллелизма" в главном гнезде циклов. | |||
Вещественный вариант обратной подстановки реализован как в основных библиотеках программ отечественных организаций, так и в западных пакетах LINPACK, LAPACK, ScaLAPACK и др. | |||
При этом в отечественных реализациях, как правило, выполнены стандартные требования к методу с точки зрения ошибок округления, то есть, реализован режим накопления, и обычно нет лишних операций. | |||
Реализация точечного варианта алгоритма в современных западных пакетах обычно происходит из одной и той же реализации метода в LINPACK, использующей библиотеку BLAS. | |||
Для большинства современных пакетов имеется также блочный вариант алгоритма обратной подстановки, в том числе и тот, граф которого топологически тождествен графу точечного варианта. Из-за того, что количество читаемых данных примерно равно количеству операций, блочность может дать некоторое ускорение работы благодаря лучшему использованию кэшей процессоров. Именно в направлении оптимизации кэширования и следует сосредоточить основные усилия при оптимизации работы программы. | |||
=== Результаты прогонов === | |||
=== Выводы для классов архитектур === | === Выводы для классов архитектур === | ||
Если исходить из структуры алгоритма, то при реализации на суперкомпьютерах следует выполнить две вещи. Во-первых, для минимизации обменов между узлами следует избрать блочный вариант, в котором или все элементы матрицы доступны на всех узлах, или заранее распределены по узлам. В таком случае количество передаваемых между узлами данных будет невелико по сравнению с количеством арифметических операций. Но при такой организации работы получится, что наибольшие временные затраты будут связаны с неоптимальностью обработки отдельных блоков. Поэтому, видимо, следует сначала оптимизировать не блочный алгоритм в целом, а подпрограммы, используемые на отдельных процессорах: точечный метод обратной подстановки, перемножения матриц и др. подпрограммы. Ниже содержится информация о возможном направлении такой оптимизации. | |||
== Литература == | |||
<references /> | |||
[[En:Backward substitution]] | |||
[[Категория:Законченные статьи]] | |||
Текущая версия от 11:33, 18 июля 2022
| Обратная подстановка для неособенной верхней треугольной матрицы | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ O(n^2) }[/math] |
| Объём входных данных | [math]\displaystyle{ O(n^2) }[/math] |
| Объём выходных данных | [math]\displaystyle{ n }[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]\displaystyle{ O(n) }[/math] |
| Ширина ярусно-параллельной формы | [math]\displaystyle{ O(n) }[/math] |
Основные авторы описания: А.В.Фролов.
Свойства и структура алгоритма
Общее описание алгоритма
Обратная подстановка - решение системы линейных алгебраических уравнений (СЛАУ) [math]\displaystyle{ Ux = y }[/math] с верхней треугольной матрицей [math]\displaystyle{ U }[/math]. Матрица [math]\displaystyle{ U }[/math] может быть одной из составляющих матрицы [math]\displaystyle{ A }[/math] в каких-либо разложениях и получается либо из [math]\displaystyle{ LU }[/math]-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, разложение Холецкого и др.), либо из других (например из QR-разложения). В силу треугольности [math]\displaystyle{ U }[/math] решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.
В[1] методом обратной подстановки назван также и метод решения СЛАУ с нижней треугольной матрицей. Там же отмечено, что в литературе иногда под обратной подстановкой имеют в виду, как и здесь, только решения СЛАУ с верхней треугольной матрицей, а решение нижних треугольных систем называют прямой подстановкой. Такой же системы названий будем придерживаться и здесь, во избежание одноимённого названия различных алгоритмов. Кроме того, обратная подстановка, представленная здесь, одновременно может быть частью метода Гаусса для решения СЛАУ, а именно - его обратным ходом, чего нельзя сказать про прямую подстановку.
Существует метод со сходным названием - Обратная подстановка с нормировкой. При том, что он решает, по существу, ту же задачу, что и простая обратная подстановка, его схема несколько сложнее. Это связано со специальными мерами по уменьшению влияния ошибок округления на результат. Алгоритм обратной подстановки с нормировкой в данном разделе не рассматривается.
Математическое описание алгоритма
Исходные данные: верхняя треугольная матрица [math]\displaystyle{ U }[/math] (элементы [math]\displaystyle{ u_{ij} }[/math]), вектор правой части [math]\displaystyle{ y }[/math] (элементы [math]\displaystyle{ y_{i} }[/math]).
Вычисляемые данные: вектор решения [math]\displaystyle{ x }[/math] (элементы [math]\displaystyle{ x_{i} }[/math]).
Формулы метода:
- [math]\displaystyle{ \begin{align} x_{n} & = y_{n}/u_{nn} \\ x_{i} & = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}, \quad i \in [1, n - 1]. \end{align} }[/math]
Существует также блочная версия метода, однако в данном описании разобран только точечный метод.
Вычислительное ядро алгоритма
Вычислительное ядро алгоритма обратной подстановки можно составить из множественных (всего их [math]\displaystyle{ n-1 }[/math]) вычислений скалярных произведений строк матрицы [math]\displaystyle{ U }[/math] (за исключением диагонального элемента) на уже вычисленную часть вектора [math]\displaystyle{ x }[/math]:
- [math]\displaystyle{ \sum_{j = i+1}^{n} u_{ij} x_{j} }[/math]
в режиме накопления или, в зависимости от требований задачи, без его использования, с их последующим вычитанием из компоненты вектора [math]\displaystyle{ y }[/math] и деления на диагональный элемент матрицы [math]\displaystyle{ U }[/math]. В отечественных реализациях, даже в последовательных версиях, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм вида
- [math]\displaystyle{ y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} }[/math]
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента вектора», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому вычислительным ядром алгоритма следует считать не вычисления скалярных произведений, а вычисления выражений
- [math]\displaystyle{ y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} }[/math]
в режиме накопления или, в зависимости от требований задачи, без его использования, плюс деления результатов этих вычислений на диагональные элементы матрицы.
Макроструктура алгоритма
Как уже упоминалось в описании ядра алгоритма, основную часть метода обратной подстановки составляют множественные (всего их [math]\displaystyle{ n-1 }[/math]) вычисления сумм
- [math]\displaystyle{ y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} }[/math]
в режиме накопления или без его использования, а также деления результатов этих вычислений на диагональные элементы матрицы.
Схема реализации последовательного алгоритма
Последовательность действий описываемого алгоритма:
1. [math]\displaystyle{ x_{n} = y_{n}/u_{nn} }[/math]
Далее для всех [math]\displaystyle{ i }[/math] от [math]\displaystyle{ n-1 }[/math] до [math]\displaystyle{ 1 }[/math] по убыванию выполняются
2. [math]\displaystyle{ x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii} }[/math].
Особо отметим, что вычисления сумм вида [math]\displaystyle{ y_{i} - \sum{}_{j = i+1}^{n} u_{ij} x_{j} }[/math] производят в режиме накопления вычитанием из [math]\displaystyle{ y_{i} }[/math] произведений [math]\displaystyle{ u_{ij} x_{j} }[/math] для [math]\displaystyle{ j }[/math] от [math]\displaystyle{ n }[/math] до [math]\displaystyle{ i + 1 }[/math], c убыванием [math]\displaystyle{ j }[/math]. Использование другого порядка выполнения суммирования приводит к резкому ухудшению параллельных свойств алгоритма, хотя, к сожалению, это кое-где встречается в литературе и пакетах программ. В качестве примера использования другого порядка вычислений можно привести фрагмент программы из[2], где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями используемого авторами книги старого диалекта Фортрана.
Последовательная сложность алгоритма
Для алгоритма обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка [math]\displaystyle{ n }[/math] в последовательном (наиболее быстром) варианте требуется:
- [math]\displaystyle{ n }[/math] делений,
- [math]\displaystyle{ \frac{n^2-n}{2} }[/math] сложений (вычитаний),
- [math]\displaystyle{ \frac{n^2-n}{2} }[/math] умножений.
Выполнение умножений и сложений (вычитаний) составляет основную часть алгоритма.
При этом использование режима накопления требует выполнения умножений и вычитаний в режиме двойной точности (или использования функции аналогичной DPROD на языке Фортран).
Таким образом, при классификации по последовательной сложности, метод обратной подстановки относится к алгоритмам со сложностью [math]\displaystyle{ O(n^2) }[/math].
Информационный граф
Опишем граф алгоритма как аналитически, так и в виде рисунка.
Граф алгоритма обратной подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности.
Первая группа вершин расположена в одномерной области, соответствующая ей операция выполняет функцию деления. Естественно введённая единственная координата каждой из вершин [math]\displaystyle{ i }[/math] изменяется в диапазоне от [math]\displaystyle{ n }[/math] до [math]\displaystyle{ 1 }[/math], принимая все целочисленные значения.
Делимое в этой операции:
- при [math]\displaystyle{ i = n }[/math] — элемент входных данных, а именно [math]\displaystyle{ y_{n} }[/math];
- при [math]\displaystyle{ i \lt n }[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]\displaystyle{ i }[/math], [math]\displaystyle{ i+1 }[/math].
Делитель для этой операции - элемент входных данных, а именно [math]\displaystyle{ u_{nn} }[/math].
Результат выполнения операции является выходным данным [math]\displaystyle{ x_{i} }[/math].
Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]\displaystyle{ a-bc }[/math]. Естественно введённые координаты области:
- [math]\displaystyle{ i }[/math] — меняется в диапазоне от [math]\displaystyle{ n-1 }[/math] до [math]\displaystyle{ 1 }[/math], принимая все целочисленные значения;
- [math]\displaystyle{ j }[/math] — меняется в диапазоне от [math]\displaystyle{ n }[/math] до [math]\displaystyle{ i+1 }[/math], принимая все целочисленные значения.
Аргументы операции:
- [math]\displaystyle{ a }[/math]:
- при [math]\displaystyle{ j = n }[/math] элемент входных данных [math]\displaystyle{ y_{i} }[/math];
- при [math]\displaystyle{ j \lt n }[/math] — результат выполнения операции, соответствующей вершине из второй группы, с координатами [math]\displaystyle{ i, j+1 }[/math];
- [math]\displaystyle{ b }[/math] — элемент входных данных, а именно [math]\displaystyle{ u_{ij} }[/math];
- [math]\displaystyle{ c }[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатой [math]\displaystyle{ j }[/math].
Результат срабатывания операции является промежуточным данным алгоритма.
Описанный граф представлен на рис. 1, выполненном для случая [math]\displaystyle{ n = 5 }[/math]. Здесь вершины первой группы обозначены жёлтым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена только подача входных данных из вектора [math]\displaystyle{ y }[/math], подача элементов матрицы [math]\displaystyle{ U }[/math], идущая во все вершины, на рисунке не представлена.
Ресурс параллелизма алгоритма
Для обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка [math]\displaystyle{ n }[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:
- [math]\displaystyle{ n }[/math] ярусов делений (в каждом из ярусов одно деление),
- по [math]\displaystyle{ n - 1 }[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n-1 }[/math].
Таким образом, в параллельном варианте, в отличие от последовательного, вычисления делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах ЯПФ отдельных делений может породить и другие проблемы. Например, при реализации метода обратной подстановки на ПЛИСах остальные вычисления (умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; деления из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени.
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.
Таким образом, при классификации по высоте ЯПФ, алгоритм обратной подстановки относится к алгоритмам с линейной сложностью. При классификации по ширине ЯПФ его сложность также будет линейной.
Входные и выходные данные алгоритма
Входные данные: верхняя треугольная матрица [math]\displaystyle{ U }[/math] (элементы [math]\displaystyle{ u_{ij} }[/math]), вектор правой части [math]\displaystyle{ y }[/math] (элементы [math]\displaystyle{ y_{i} }[/math]).
Объём входных данных: [math]\displaystyle{ \frac{n (n + 3)}{2} }[/math] (в силу треугольности достаточно хранить только ненулевые элементы матрицы [math]\displaystyle{ U }[/math]).
Выходные данные: вектор решения [math]\displaystyle{ x }[/math] (элементы [math]\displaystyle{ x_{i} }[/math]).
Объём выходных данных: [math]\displaystyle{ n }[/math].
Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейным (отношение квадратической к линейной).
При этом вычислительная мощность алгоритма обратной подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь константа.
Представленный алгоритм обратной подстановки является полностью детерминированным. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и меняет сложность с линейной на квадратичную.
Наличие линейного количества ярусов ЯПФ, состоящих из одного-единственного деления, потенциально замедляющее параллельные реализации алгоритма, является его характерным "узким местом", особенно в сравнении со схожей по решаемой математической задаче прямой подстановке, где диагональные элементы единичны. В связи с этим для решения СЛАУ предпочтительны такие разложения, содержащие треугольные матрицы, где в треугольных матрицах диагональные элементы единичны. В тех же случаях, когда получаются неособенные треугольные матрицы, их желательно предварительно, до решения СЛАУ с ними, преобразовать в произведение диагональной и треугольной с единичными диагональными элементами.
У алгоритма обратной подстановки существует несколько блочных вариантов. Граф некоторых из них совпадает с графом точечного варианта, различия связаны в основном с порядком прохождения основных циклов алгоритма, а именно, с их развёртыванием и перестановкой. Эти приёмы могут помочь в оптимизации обменов на конкретных вычислительных системах.
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
В простейшем варианте метод обратной подстановки на Фортране можно записать так:
X(N) = Y(N) / U(N,N)
DO I = N-1, 1, -1
S = Y(I)
DO J = N, I+1, -1
S = S - DPROD(U(I,J), X(J))
END DO
X(I) = S / U(I,I)
END DO
При этом для реализации режима накопления переменная [math]\displaystyle{ S }[/math] должна быть двойной точности.
Возможные способы и особенности параллельной реализации алгоритма
Вариантов параллельной реализации алгоритма не так уж и много, если не использовать то, что оба главных цикла можно развернуть, перейдя, таким образом, к блочной версии. Версии без развёртывания циклов возможны как с полностью параллельными циклами по I:
DO PARALLEL I = 1, N
X(I) = Y(I)
END DO
DO J = N, 1, -1
X(J) = X(J) / U(J,J)
DO PARALLEL I = 1, J-1
X(I) = X(I) - U(I,J)*X(J)
END DO
END DO
так и с использованием "скошенного параллелизма" в главном гнезде циклов.
Вещественный вариант обратной подстановки реализован как в основных библиотеках программ отечественных организаций, так и в западных пакетах LINPACK, LAPACK, ScaLAPACK и др. При этом в отечественных реализациях, как правило, выполнены стандартные требования к методу с точки зрения ошибок округления, то есть, реализован режим накопления, и обычно нет лишних операций. Реализация точечного варианта алгоритма в современных западных пакетах обычно происходит из одной и той же реализации метода в LINPACK, использующей библиотеку BLAS.
Для большинства современных пакетов имеется также блочный вариант алгоритма обратной подстановки, в том числе и тот, граф которого топологически тождествен графу точечного варианта. Из-за того, что количество читаемых данных примерно равно количеству операций, блочность может дать некоторое ускорение работы благодаря лучшему использованию кэшей процессоров. Именно в направлении оптимизации кэширования и следует сосредоточить основные усилия при оптимизации работы программы.
Результаты прогонов
Выводы для классов архитектур
Если исходить из структуры алгоритма, то при реализации на суперкомпьютерах следует выполнить две вещи. Во-первых, для минимизации обменов между узлами следует избрать блочный вариант, в котором или все элементы матрицы доступны на всех узлах, или заранее распределены по узлам. В таком случае количество передаваемых между узлами данных будет невелико по сравнению с количеством арифметических операций. Но при такой организации работы получится, что наибольшие временные затраты будут связаны с неоптимальностью обработки отдельных блоков. Поэтому, видимо, следует сначала оптимизировать не блочный алгоритм в целом, а подпрограммы, используемые на отдельных процессорах: точечный метод обратной подстановки, перемножения матриц и др. подпрограммы. Ниже содержится информация о возможном направлении такой оптимизации.