Прямая подстановка (вещественный вариант)
Описание свойств и структуры алгоритма
Словесное описание алгоритма
Прямая подстановка - решение СЛАУ [math]\displaystyle{ Lx = y }[/math] с левой треугольной матрицей [math]\displaystyle{ L }[/math]. Матрица [math]\displaystyle{ L }[/math] - одна из составляющих матрицы [math]\displaystyle{ A }[/math] и получается либо из [math]\displaystyle{ LU }[/math]-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса), либо из других разложений. В силу треугольности [math]\displaystyle{ L }[/math] решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.
В [1] метод решения СЛАУ с левой треугольной матрицей назван методом обратной подстановки. Там же отмечено, что в литературе иногда под обратной подстановкой имеют в виду, как и здесь, только решения СЛАУ с правой треугольной матрицей, а решение левых треугольных систем называют прямой подстановкой. Такой же системы названий будем придерживаться и мы, во избежание одноимённого названия разных алгоритмов. Кроме того, обратная подстановка, представленная в этой энциклопедии алгоритмов, одновременно является частью метода Гаусса для решения СЛАУ, а именно - его обратным ходом, чего нельзя сказать про представленную здесь прямую подстановку.
Общая структура прямой подстановки с неособенной левой треугольной матрицей, тем не менее, практически полностью совпадает со структурой обратной подстановки. Поэтому здесь мы рассмотрим случай, когда матрица [math]\displaystyle{ L }[/math], как полученная из разложения типа Гаусса, имеет единичные диагональные элементы.
Математическое описание
Исходные данные: левая треугольная матрица [math]\displaystyle{ L }[/math] (элементы [math]\displaystyle{ l_{ij} }[/math]), вектор правой части [math]\displaystyle{ b }[/math] (элементы [math]\displaystyle{ b_{i} }[/math]).
Вычисляемые данные: вектор решения [math]\displaystyle{ y }[/math] (элементы [math]\displaystyle{ y_{i} }[/math]).
Формулы метода:
- [math]\displaystyle{ \begin{align} y_{1} & = b_{1} \\ y_{i} & = b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j}, \quad i \in [2, n]. \end{align} }[/math]
Существует также блочная версия метода, однако в данном описании разобран только точечный метод.
Вычислительное ядро алгоритма
Вычислительное ядро прямой подстановки можно составить из множественных (всего их [math]\displaystyle{ n-1 }[/math]) вычислений скалярных произведений строк матрицы [math]\displaystyle{ L }[/math] на уже вычисленную часть вектора [math]\displaystyle{ y }[/math]:
- [math]\displaystyle{ \sum_{j = 1}^{i-1} l_{ij} y_{j} }[/math]
в режиме накопления или без него, в зависимости от требований задачи, с их последующим вычитанием из компоненты вектора [math]\displaystyle{ b }[/math] и деления на диагональный элемент матрицы [math]\displaystyle{ L }[/math]. В отечественных реализациях, даже в последовательных, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм типа
- [math]\displaystyle{ b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} }[/math]
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода прямой подстановки не вычисления скалярных произведений, а вычисления выражений
- [math]\displaystyle{ b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} }[/math]
в режиме накопления или без него, в зависимости от требований задачи.
Макроструктура алгоритма
Как уже записано в описании ядра алгоритма, основную часть метода прямой подстановки составляют множественные (всего [math]\displaystyle{ n-1 }[/math]) вычисления сумм
- [math]\displaystyle{ b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} }[/math]
в режиме накопления или без него.
Описание схемы реализации последовательного алгоритма
Чтобы понять последовательность исполнения, перепишем формулы метода так:
1. [math]\displaystyle{ y_{1} = b_{1} }[/math]
Далее для всех [math]\displaystyle{ i }[/math] от [math]\displaystyle{ 2 }[/math] до [math]\displaystyle{ n }[/math] по возрастанию выполняются
2. [math]\displaystyle{ y_{i} = b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} }[/math]
Особо отметим, что вычисления сумм вида [math]\displaystyle{ b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} }[/math] производят в режиме накопления вычитанием из [math]\displaystyle{ b_{i} }[/math] произведений [math]\displaystyle{ l_{ij} y_{j} }[/math] для [math]\displaystyle{ j }[/math] от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ i-1 }[/math], c возрастанием [math]\displaystyle{ j }[/math]. Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма.
Последовательная сложность алгоритма
Для прямой подстановки порядка n в последовательном (наиболее быстром) варианте требуется:
- по [math]\displaystyle{ \frac{n^2-n}{2} }[/math] умножений и сложений (вычитаний).
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает затраты во времени, требуемом для выполнения прямой подстановки.
При классификации по последовательной сложности, таким образом, метод обратной подстановки относится к алгоритмам с квадратической сложностью.
Информационный граф
Опишем граф алгоритма как аналитически, так и в виде рисунка.
Граф алгоритма прямой подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности.
Описание ресурса параллелизма алгоритма
Для прямой подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:
- по [math]\displaystyle{ n - 1 }[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ n-1 }[/math].
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.
При классификации по высоте ЯПФ, таким образом, метод прямой подстановки относится к алгоритмам с линейной сложностью. При классификации по ширине ЯПФ его сложность также будет линейной.
Описание входных и выходных данных
Входные данные: левая треугольная матрица [math]\displaystyle{ L }[/math] (элементы [math]\displaystyle{ l_{ij} }[/math]), вектор правой части [math]\displaystyle{ b }[/math] (элементы [math]\displaystyle{ b_{i} }[/math]).
Объём входных данных: [math]\displaystyle{ \frac{n (n + 1)}{2} }[/math] (в силу треугольности и единичности диагональных элементов достаточно хранить только поддиагональные элементы матрицы [math]\displaystyle{ L }[/math]).
Выходные данные: вектор решения [math]\displaystyle{ y }[/math] (элементы [math]\displaystyle{ y_{i} }[/math]).
Объём выходных данных: [math]\displaystyle{ n }[/math].
Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейным (отношение квадратической к линейной).
При этом вычислительная мощность алгоритма прямой подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь константа.
При этом алгоритм прямой подстановки полностью детерминирован. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и делает параллельную сложность квадратической.
Программная реализация
Особенности реализации последовательного алгоритма
В простейшем варианте метод прямой подстановки на Фортране можно записать так:
Y(1) = B(1)
DO I = 2, N-1
S = B(I)
DO J = 1, I-1
S = S - DPROD(L(I,J), Y(J))
END DO
END DO
При этом для реализации режима накопления переменная [math]\displaystyle{ S }[/math] должна быть двойной точности.
Литература
- В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. М.: Наука, 1984, стр. 182.