Обратная подстановка (вещественный вариант): различия между версиями
[непроверенная версия] | [непроверенная версия] |
Frolov (обсуждение | вклад) |
Frolov (обсуждение | вклад) |
||
Строка 57: | Строка 57: | ||
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} - \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>. '''''Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма''''', хотя, к сожалению, остаются кое-где в литературе и пакетах программ. | + | Особо отметим, что вычисления сумм вида <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>. '''''Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма''''', хотя, к сожалению, остаются кое-где в литературе и пакетах программ.В качестве примера такого порядка можно привести фрагмент программы из [2], где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями старого диалекта Фортрана. |
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === |
Версия 15:52, 17 сентября 2014
Содержание
- 1 Описание свойств и структуры алгоритма
- 1.1 Словесное описание алгоритма
- 1.2 Математическое описание
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Описание схемы реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Описание ресурса параллелизма алгоритма
- 1.9 Описание входных и выходных данных
- 1.10 Свойства алгоритма
- 2 Программная реализация
- 3 Литература
1 Описание свойств и структуры алгоритма
1.1 Словесное описание алгоритма
Обратный ход метода Гаусса, или обратная подстановка - решение СЛАУ [math]Ux = y[/math] с правой треугольной матрицей [math]U[/math]. Матрица [math]U[/math] - одна из составляющих матрицы [math]A[/math] и получается либо из [math]LU[/math]-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, разложение Холецкого и др.), либо из других разложений. В силу треугольности [math]U[/math] решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.
В [1] методом обратной подстановки назван также и метод решения СЛАУ с левой треугольной матрицей. Там же отмечено, что в литературе иногда под обратной подстановкой имеют в виду, как и здесь, только решения СЛАУ с правой треугольной матрицей, а решение левых треугольных систем называют прямой подстановкой. Такой же системы названий будем придерживаться и мы, во избежание одноимённого названия разных алгоритмов. Кроме того, обратная подстановка, представленная здесь, одновременно является частью метода Гаусса для решения СЛАУ, а именно - его обратным ходом, чего нельзя сказать про прямую подстановку.
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]. Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма, хотя, к сожалению, остаются кое-где в литературе и пакетах программ.В качестве примера такого порядка можно привести фрагмент программы из [2], где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями старого диалекта Фортрана.
1.6 Последовательная сложность алгоритма
Для обратной подстановки порядка n в последовательном (наиболее быстром) варианте требуется:
- [math]n[/math] делений,
- по [math]\frac{n^2-n}{2}[/math] умножений и сложений (вычитаний) — основная часть алгоритма.
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения обратного хода метода Гаусса.
При классификации по последовательной сложности, таким образом, метод обратной подстановки относится к алгоритмам с квадратической сложностью.
1.7 Информационный граф
Опишем граф алгоритма как аналитически, так и в виде рисунка.
Граф алгоритма обратной подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности.
Первая группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию деления. Естественно введённая единственная координата каждой из вершин [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];
Результат срабатывания операции является промежуточным данным алгоритма.
Описанный граф можно посмотреть на рисунке, выполненном для случая [math]n = 5[/math]. Здесь вершины первой группы обозначены розовым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена подача только входных данных из вектора [math]y_{i}[/math], подача элементов матрицы [math]U[/math], идущая во все вершины, на рисунке не представлена.
1.8 Описание ресурса параллелизма алгоритма
Для обратной подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:
- [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
END DO
При этом для реализации режима накопления переменная [math]S[/math] должна быть двойной точности.
3 Литература
- В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. - М.: Наука, 1984.
- Дж.Форсайт, К.Моллер. Численное решение систем линейных алгебраических уравнений. - М.:Мир, 1969.