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

Прямая подстановка для нижней треугольной матрицы с единичной диагональю
Последовательный алгоритм
Последовательная сложность	[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{ Lx = y }[/math] с нижней треугольной матрицей [math]\displaystyle{ L }[/math]. Матрица [math]\displaystyle{ L }[/math] - одна из составляющих матрицы [math]\displaystyle{ A }[/math] и получается либо из [math]\displaystyle{ LU }[/math]-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса), либо из других разложений. В силу треугольности [math]\displaystyle{ L }[/math] решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.

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

Общая структура прямой подстановки с неособенной нижней треугольной матрицей, тем не менее, практически полностью совпадает со структурой обратной подстановки. Поэтому здесь мы рассмотрим случай, когда матрица [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 в Фортране), что ещё больше увеличивает затраты во времени, требуемом для выполнения прямой подстановки.

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

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

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

Рисунок 1. Прямая подстановка

Граф алгоритма прямой подстановки состоит из одной группы вершин, расположенной в целочисленных узлах двумерной области, соответствующая ей операция [math]\displaystyle{ a-bc }[/math].

Естественно введённые координаты области таковы:

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

Аргументы операции следующие:

• [math]\displaystyle{ a }[/math]:
  • при [math]\displaystyle{ j = 1 }[/math] элемент входных данных [math]\displaystyle{ b_{i} }[/math];
  • при [math]\displaystyle{ j \gt 1 }[/math] — результат срабатывания операции, соответствующей вершине с координатами [math]\displaystyle{ i, j-1 }[/math];
• [math]\displaystyle{ b }[/math] — элемент входных данных, а именно [math]\displaystyle{ l_{ij} }[/math];
• [math]\displaystyle{ c }[/math]:
  • при [math]\displaystyle{ j = 1 }[/math] элемент входных данных [math]\displaystyle{ b_{1} }[/math];
  • при [math]\displaystyle{ j \gt 1 }[/math] — результат срабатывания операции, соответствующей вершине с координатой [math]\displaystyle{ j, j-1 }[/math];

Результат срабатывания операции является:

• при [math]\displaystyle{ j \lt i-1 }[/math] - промежуточным данным алгоритма;
• при [math]\displaystyle{ j = i-1 }[/math] - выходным данным.

Описанный граф можно посмотреть на рисунке, выполненном для случая [math]\displaystyle{ n = 5 }[/math]. Здесь вершины обозначены зелёным цветом. Подача входных данных из вектора [math]\displaystyle{ b }[/math], идущая в вершины левого столбца, и матрицы [math]\displaystyle{ L }[/math], идущая во все вершины, на рисунке не представлена.

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

Для прямой подстановки порядка 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] должна быть двойной точности.