Обратная подстановка (вещественный вариант)

Обратная подстановка для неособенной верхней треугольной матрицы
Последовательный алгоритм
Последовательная сложность	[math]O(n^2)[/math]
Объём входных данных	[math]O(n^2)[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(n)[/math]
Ширина ярусно-параллельной формы	[math]O(n)[/math]

Основные авторы описания: А.В.Фролов.

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Обратная подстановка - решение системы линейных алгебраических уравнений (СЛАУ) [math]Ux = y[/math] с верхней треугольной матрицей [math]U[/math]. Матрица [math]U[/math] может быть одной из составляющих матрицы [math]A[/math] в каких-либо разложениях и получается либо из [math]LU[/math]-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, разложение Холецкого и др.), либо из других (например из QR-разложения). В силу треугольности [math]U[/math] решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.

В методом обратной подстановки назван также и метод решения СЛАУ с нижней треугольной матрицей. Там же отмечено, что в литературе иногда под обратной подстановкой имеют в виду, как и здесь, только решения СЛАУ с верхней треугольной матрицей, а решение нижних треугольных систем называют прямой подстановкой. Такой же системы названий будем придерживаться и здесь, во избежание одноимённого названия различных алгоритмов. Кроме того, обратная подстановка, представленная здесь, одновременно может быть частью метода Гаусса для решения СЛАУ, а именно - его обратным ходом, чего нельзя сказать про прямую подстановку.

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

1.2 Математическое описание алгоритма

Исходные данные: верхняя треугольная матрица [math]U[/math] (элементы [math]u_{ij}[/math]), вектор правой части [math]y[/math] (элементы [math]y_{i}[/math]).

Вычисляемые данные: вектор решения [math]x[/math] (элементы [math]x_{i}[/math]).

Формулы метода:

[math] \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]

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

1.3 Вычислительное ядро алгоритма

Вычислительное ядро алгоритма обратной подстановки можно составить из множественных (всего их [math]n-1[/math]) вычислений скалярных произведений строк матрицы [math]U[/math] (за исключением диагонального элемента) на уже вычисленную часть вектора [math]x[/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]

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

1.4 Макроструктура алгоритма

Как уже упоминалось в описании ядра алгоритма, основную часть метода обратной подстановки составляют множественные (всего их [math]n-1[/math]) вычисления сумм

[math]y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

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

1.5 Схема реализации последовательного алгоритма

Последовательность действий описываемого алгоритма:

1. [math]x_{n} = y_{n}/u_{nn}[/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].

Особо отметим, что вычисления сумм вида [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]. Использование другого порядка выполнения суммирования приводит к резкому ухудшению параллельных свойств алгоритма, хотя, к сожалению, это кое-где встречается в литературе и пакетах программ. В качестве примера использования другого порядка вычислений можно привести фрагмент программы из, где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями используемого авторами книги старого диалекта Фортрана.

1.6 Последовательная сложность алгоритма

Для алгоритма обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка [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].

1.7 Информационный граф

Опишем граф алгоритма как аналитически, так и в виде рисунка.

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

Рис. 1. Обратная подстановка

Первая группа вершин расположена в одномерной области, соответствующая ей операция выполняет функцию деления. Естественно введённая единственная координата каждой из вершин [math]i[/math] изменяется в диапазоне от [math]n[/math] до [math]1[/math], принимая все целочисленные значения.

Делимое в этой операции:

• при [math]i = n[/math] — элемент входных данных, а именно [math]y_{n}[/math];
• при [math]i \lt n[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]i[/math], [math]i+1[/math].

Делитель для этой операции - элемент входных данных, а именно [math]u_{nn}[/math].

Результат выполнения операции является выходным данным [math]x_{i}[/math].

Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]a-bc[/math]. Естественно введённые координаты области:

• [math]i[/math] — меняется в диапазоне от [math]n-1[/math] до [math]1[/math], принимая все целочисленные значения;
• [math]j[/math] — меняется в диапазоне от [math]n[/math] до [math]i+1[/math], принимая все целочисленные значения.

Аргументы операции:

• [math]a[/math]:
  • при [math]j = n[/math] элемент входных данных [math]y_{i}[/math];
  • при [math]j \lt n[/math] — результат выполнения операции, соответствующей вершине из второй группы, с координатами [math]i, j+1[/math];
• [math]b[/math] — элемент входных данных, а именно [math]u_{ij}[/math];
• [math]c[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатой [math]j[/math].

Результат срабатывания операции является промежуточным данным алгоритма.

Описанный граф представлен на рис. 1, выполненном для случая [math]n = 5[/math]. Здесь вершины первой группы обозначены жёлтым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена только подача входных данных из вектора [math]y[/math], подача элементов матрицы [math]U[/math], идущая во все вершины, на рисунке не представлена.

1.8 Ресурс параллелизма алгоритма

Для обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка [math]n[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

• [math]n[/math] ярусов делений (в каждом из ярусов одно деление),
• по [math]n - 1[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от [math]1[/math] до [math]n-1[/math].

Таким образом, в параллельном варианте, в отличие от последовательного, вычисления делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах ЯПФ отдельных делений может породить и другие проблемы. Например, при реализации метода обратной подстановки на ПЛИСах остальные вычисления (умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; деления из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени.

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

Таким образом, при классификации по высоте ЯПФ, алгоритм обратной подстановки относится к алгоритмам с линейной сложностью. При классификации по ширине ЯПФ его сложность также будет линейной.

1.9 Входные и выходные данные алгоритма

Входные данные: верхняя треугольная матрица [math]U[/math] (элементы [math]u_{ij}[/math]), вектор правой части [math]y[/math] (элементы [math]y_{i}[/math]).

Объём входных данных: [math]\frac{n (n + 3)}{2}[/math] (в силу треугольности достаточно хранить только ненулевые элементы матрицы [math]U[/math]).

Выходные данные: вектор решения [math]x[/math] (элементы [math]x_{i}[/math]).

Объём выходных данных: [math]n[/math].

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейным (отношение квадратической к линейной).

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

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

Наличие линейного количества ярусов ЯПФ, состоящих из одного-единственного деления, потенциально замедляющее параллельные реализации алгоритма, является его характерным "узким местом", особенно в сравнении со схожей по решаемой математической задаче прямой подстановке, где диагональные элементы единичны. В связи с этим для решения СЛАУ предпочтительны такие разложения, содержащие треугольные матрицы, где в треугольных матрицах диагональные элементы единичны. В тех же случаях, когда получаются неособенные треугольные матрицы, их желательно предварительно, до решения СЛАУ с ними, преобразовать в произведение диагональной и треугольной с единичными диагональными элементами.

У алгоритма обратной подстановки существует несколько блочных вариантов. Граф некоторых из них совпадает с графом точечного варианта, различия связаны в основном с порядком прохождения основных циклов алгоритма, а именно, с их развёртыванием и перестановкой. Эти приёмы могут помочь в оптимизации обменов на конкретных вычислительных системах.

2 Программная реализация алгоритма

2.1 Особенности реализации последовательного алгоритма

В простейшем варианте метод обратной подстановки на Фортране можно записать так:

        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]S[/math] должна быть двойной точности.

2.2 Возможные способы и особенности параллельной реализации алгоритма

Вариантов параллельной реализации алгоритма не так уж и много, если не использовать то, что оба главных цикла можно развернуть, перейдя, таким образом, к блочной версии. Версии без развёртывания циклов возможны как с полностью параллельными циклами по 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.

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

2.4 Выводы для классов архитектур

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