Уровень алгоритма

Прямая подстановка (вещественный вариант): различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][досмотренная версия]
(Новая страница: «== Описание свойств и структуры алгоритма == === Словесное описание алгоритма === '''Прямая п…»)
 
 
(не показано 27 промежуточных версий 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>Lx = y</math> с левой треугольной матрицей <math>L</math>. Матрица <math>L</math> - одна из составляющих матрицы <math>A</math> и получается либо из <math>LU</math>-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса), либо из других разложений. В силу треугольности <math>L</math> решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.
+
== Свойства и структура алгоритма ==
  
В [1] метод решения СЛАУ с ''левой треугольной матрицей'' назван методом '''обратной подстановки'''. Там же отмечено, что в литературе иногда под ''обратной подстановкой'' имеют в виду, как и здесь, только решения СЛАУ с ''правой треугольной матрицей'', а решение ''левых'' треугольных систем называют прямой подстановкой. Такой же системы названий будем придерживаться и мы, во избежание одноимённого названия разных алгоритмов. Кроме того, [[Обратная подстановка (вещественный вариант)|обратная подстановка]], представленная в этой энциклопедии алгоритмов, одновременно является частью '''метода Гаусса для решения СЛАУ''', а именно - его '''обратным ходом''', чего нельзя сказать про представленную здесь прямую подстановку.
+
=== Общее описание алгоритма ===
  
Общая структура прямой подстановки с неособенной левой треугольной матрицей, тем не менее, практически полностью совпадает со структурой [[Обратная подстановка (вещественный вариант)|обратной подстановки]]. Поэтому здесь мы рассмотрим случай, когда матрица <math>L</math>, как полученная из разложения типа Гаусса, имеет единичные диагональные элементы.  
+
'''Прямая подстановка''' - решение СЛАУ <math>Lx = y</math> с нижней треугольной матрицей <math>L</math>. Матрица <math>L</math> - одна из составляющих матрицы <math>A</math> и получается либо из <math>LU</math>-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса), либо из других разложений. В силу треугольности <math>L</math> решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.
  
=== Математическое описание ===
+
В<ref>В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. М.: Наука, 1984, стр. 182.</ref> метод решения СЛАУ с ''левой треугольной матрицей'' назван методом '''обратной подстановки'''. Там же отмечено, что в литературе иногда под ''обратной подстановкой'' имеют в виду, как и здесь, только решения СЛАУ с ''правой треугольной матрицей'', а решение ''левых'' треугольных систем называют прямой подстановкой. Такой же системы названий будем придерживаться и мы, во избежание одноимённого названия разных алгоритмов. Кроме того, [[Обратная подстановка (вещественный вариант)|обратная подстановка]], представленная в этой энциклопедии алгоритмов, одновременно является частью '''метода Гаусса для решения СЛАУ''', а именно - его '''обратным ходом''', чего нельзя сказать про представленную здесь прямую подстановку.
  
Исходные данные: левая треугольная матрица <math>L</math> (элементы <math>l_{ij}</math>), вектор правой части <math>b</math> (элементы <math>b_{i}</math>).
+
Общая структура прямой подстановки с неособенной нижней треугольной матрицей, тем не менее, практически полностью совпадает со структурой [[Обратная подстановка (вещественный вариант)|обратной подстановки]]. Поэтому здесь мы рассмотрим случай, когда матрица <math>L</math>, как полученная из разложения типа Гаусса, имеет единичные диагональные элементы.
 +
 
 +
=== Математическое описание алгоритма ===
 +
 
 +
Исходные данные: нижняя треугольная матрица <math>L</math> (элементы <math>l_{ij}</math>), вектор правой части <math>b</math> (элементы <math>b_{i}</math>).
  
 
Вычисляемые данные: вектор решения <math>y</math> (элементы <math>y_{i}</math>).
 
Вычисляемые данные: вектор решения <math>y</math> (элементы <math>y_{i}</math>).
Строка 27: Строка 38:
 
=== Вычислительное ядро алгоритма ===
 
=== Вычислительное ядро алгоритма ===
  
Вычислительное ядро обратной подстановки можно составить из множественных (всего их <math>n-1</math>) вычислений скалярных произведений строк матрицы <math>L</math> на уже вычисленную часть вектора <math>y</math>:
+
Вычислительное ядро прямой подстановки можно составить из множественных (всего их <math>n-1</math>) вычислений скалярных произведений строк матрицы <math>L</math> на уже вычисленную часть вектора <math>y</math>:
  
 
:<math> \sum_{j = 1}^{i-1} l_{ij} y_{j} </math>  
 
:<math> \sum_{j = 1}^{i-1} l_{ij} y_{j} </math>  
Строка 35: Строка 46:
 
:<math> b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} </math>
 
:<math> b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} </math>
  
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода не вычисления скалярных произведений, а вычисления выражений
+
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода прямой подстановки не вычисления скалярных произведений, а вычисления выражений
  
 
:<math> b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} </math>
 
:<math> b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} </math>
Строка 43: Строка 54:
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
  
Как уже записано в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть метода обратной подстановки составляют множественные (всего <math>n-1</math>) вычисления сумм
+
Как уже записано в [[#Вычислительное ядро алгоритма|описании ядра алгоритма]], основную часть метода прямой подстановки составляют множественные (всего <math>n-1</math>) вычисления сумм
  
:<math>y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} </math>
+
:<math> b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} </math>
  
в режиме накопления или без него, плюс деления результатов этих вычислений на диагональные элементы матрицы.
+
в режиме накопления или без него.
  
=== Описание схемы реализации последовательного алгоритма ===
+
=== Схема реализации последовательного алгоритма ===
  
 
Чтобы понять последовательность исполнения, перепишем формулы метода так:
 
Чтобы понять последовательность исполнения, перепишем формулы метода так:
  
1. <math>x_{n} = y_{n}/u_{nn}</math>
+
1. <math>y_{1} = b_{1}</math>
  
Далее для всех <math>i</math> от <math>n-1</math> до <math>1</math> по убыванию выполняются
+
Далее для всех <math>i</math> от <math>2</math> до <math>n</math> по возрастанию выполняются
  
2. <math>x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}</math>
+
2. <math>y_{i} = b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{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>. '''''Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма''''', хотя, к сожалению, остаются кое-где в литературе и пакетах программ.
+
Особо отметим, что вычисления сумм вида <math>b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j}</math> производят в режиме накопления вычитанием из <math>b_{i}</math> произведений <math>l_{ij} y_{j}</math> для <math>j</math> от <math>1</math> до <math>i-1</math>, '''c возрастанием''' <math>j</math>. '''''Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма'''''.
  
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
  
Для обратной подстановки порядка n  в последовательном (наиболее быстром) варианте требуется:
+
Для прямой подстановки порядка n  в последовательном (наиболее быстром) варианте требуется:
 
 
* <math>n</math> делений,
+
* по <math>\frac{n^2-n}{2}</math> умножений и сложений (вычитаний).
* по <math>\frac{n^2-n}{2}</math> умножений и сложений (вычитаний) — ''основная часть алгоритма''.
 
 
   
 
   
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения обратного хода метода Гаусса.
+
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает затраты во времени, требуемом для выполнения прямой подстановки.
  
 
При классификации по последовательной сложности, таким образом, метод обратной подстановки относится к алгоритмам ''с квадратической сложностью''.
 
При классификации по последовательной сложности, таким образом, метод обратной подстановки относится к алгоритмам ''с квадратической сложностью''.
Строка 76: Строка 86:
 
Опишем [[глоссарий#Граф алгоритма|граф алгоритма]] как аналитически, так и в виде рисунка.  
 
Опишем [[глоссарий#Граф алгоритма|граф алгоритма]] как аналитически, так и в виде рисунка.  
  
Граф алгоритма обратной подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности.
+
[[Файл:DirectL.png|450px|thumb|left|Рисунок 1. Прямая подстановка]]
 
 
[[Файл:DirectU.png|450px|thumb|left|описание]]
 
 
 
'''Первая''' группа вершин расположена в одномерной области, соответствующая ей операция вычисляет функцию деления.
 
Естественно введённая единственная координата каждой из вершин <math>i</math> меняется в диапазоне от <math>n</math> до <math>1</math>, принимая все целочисленные значения.
 
 
 
Делимое в этой операции:
 
 
* при <math>i = n</math> — элемент ''входных данных'', а именно  <math>y_{n}</math>;
 
* при <math>i < n</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>i</math>, <math>i+1</math>.
 
 
 
Делитель для этой операции - элемент ''входных данных'', а именно  <math>u_{nn}</math>.
 
  
Результат срабатывания операции является ''выходным данным'' <math>x_{i}</math>.
+
Граф алгоритма прямой подстановки состоит из одной группы вершин, расположенной в целочисленных узлах двумерной области, соответствующая ей операция  <math>a-bc</math>.  
  
'''Вторая''' группа вершин расположена в двумерной области, соответствующая ей операция  <math>a-bc</math>.
 
 
Естественно введённые координаты области таковы:  
 
Естественно введённые координаты области таковы:  
* <math>i</math> — меняется в диапазоне от <math>n-1</math> до <math>1</math>, принимая все целочисленные значения;
+
* <math>i</math> — меняется в диапазоне от <math>2</math> до <math>n</math>, принимая все целочисленные значения;
* <math>j</math> — меняется в диапазоне от <math>n</math> до <math>i+1</math>, принимая все целочисленные значения.
+
* <math>j</math> — меняется в диапазоне от <math>1</math> до <math>i-1</math>, принимая все целочисленные значения.
  
 
Аргументы операции следующие:
 
Аргументы операции следующие:
 
*<math>a</math>:
 
*<math>a</math>:
** при <math>j = n</math> элемент входных данных <math>y_{i}</math>;  
+
** при <math>j = 1</math> элемент входных данных <math>b_{i}</math>;  
** при <math>j < n</math> — результат срабатывания операции, соответствующей вершине из второй группы, с координатами <math>i, j+1</math>;
+
** при <math>j > 1</math> — результат срабатывания операции, соответствующей вершине с координатами <math>i, j-1</math>;
*<math>b</math> — элемент ''входных данных'', а именно  <math>u_{ij}</math>;
+
*<math>b</math> — элемент ''входных данных'', а именно  <math>l_{ij}</math>;
*<math>c</math> — результат срабатывания операции, соответствующей вершине из первой группы, с координатой <math>j</math>;
+
*<math>c</math>:
 +
** при <math>j = 1</math> элемент входных данных <math>b_{1}</math>;
 +
** при <math>j > 1</math> — результат срабатывания операции, соответствующей вершине с координатой <math>j, j-1</math>;
 +
 
 +
Результат срабатывания операции является:
 +
* при <math>j < i-1</math> - ''промежуточным данным'' алгоритма;
 +
* при <math>j = i-1</math> - выходным данным.
  
Результат срабатывания операции является ''промежуточным данным'' алгоритма.
+
Описанный граф можно посмотреть на рисунке, выполненном для случая <math>n = 5</math>. Здесь вершины обозначены зелёным цветом. Подача входных данных из вектора <math>b</math>, идущая в вершины левого столбца, и матрицы <math>L</math>, идущая во все вершины, на рисунке не представлена.
  
Описанный граф можно посмотреть на рисунке, выполненном для случая <math>n = 5</math>. Здесь вершины первой группы обозначены розовым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена подача только входных данных из вектора <math>y_{i}</math>, подача элементов матрицы <math>U</math>, идущая во все вершины, на рисунке не представлена.
+
=== Ресурс параллелизма алгоритма ===
  
=== Описание ресурса параллелизма алгоритма ===
+
Для прямой подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:
  
Для обратной подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:
 
* <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>.
 
   
 
   
Таким образом, в параллельном варианте, в отличие от последовательного, вычисления делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах [[глоссарий#Ярусно-параллельная форма графа алгоритма|ЯПФ]] отдельных делений может породить и другие проблемы. Например, при реализации метода обратной подстановки на ПЛИСах остальные вычисления (умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; деления из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени.
 
 
 
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.
 
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.
  
При классификации по высоте ЯПФ, таким образом, метод обратной подстановки относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''линейной''.
+
При классификации по высоте ЯПФ, таким образом, метод прямой подстановки относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''линейной''.
  
=== Описание входных и выходных данных ===
+
=== Входные и выходные данные алгоритма ===
  
'''Входные данные''': верхняя треугольная матрица <math>U</math> (элементы <math>u_{ij}</math>), вектор правой части <math>y</math> (элементы <math>y_{i}</math>).
+
'''Входные данные''': нижняя треугольная матрица <math>L</math> (элементы <math>l_{ij}</math>), вектор правой части <math>b</math> (элементы <math>b_{i}</math>).
  
'''Объём входных данных''': <math>\frac{n (n + 3)}{2}</math> (в силу треугольности достаточно хранить только диагональ и наддиагональные элементы матрицы <math>U</math>).  
+
'''Объём входных данных''': <math>\frac{n (n + 1)}{2}</math> (в силу треугольности и единичности диагональных элементов достаточно хранить только поддиагональные элементы матрицы <math>L</math>).  
  
'''Выходные данные''': вектор решения <math>x</math> (элементы <math>x_{i}</math>).
+
'''Выходные данные''': вектор решения <math>y</math> (элементы <math>y_{i}</math>).
  
'''Объём выходных данных''': <math>n</math>.
+
'''Объём выходных данных''': <math>n~.</math>
  
 
=== Свойства алгоритма ===
 
=== Свойства алгоритма ===
Строка 134: Строка 133:
 
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является ''линейным'' (отношение квадратической к линейной).  
 
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является ''линейным'' (отношение квадратической к линейной).  
  
При этом вычислительная мощность алгоритма обратной подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь ''константа''.
+
При этом вычислительная мощность алгоритма прямой подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь ''константа''.
  
При этом алгоритм обратной подстановки полностью детерминирован. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и делает параллельную сложность квадратической.
+
При этом алгоритм прямой подстановки полностью детерминирован. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и делает параллельную сложность квадратической.
  
== Программная реализация ==
+
== Программная реализация алгоритма ==
  
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
  
В простейшем варианте метод обратной подстановки на Фортране можно записать так:
+
В простейшем варианте метод прямой подстановки на Фортране можно записать так:
  
 
<source lang="fortran">
 
<source lang="fortran">
         X(N) = Y(N) / U (N, N)
+
         Y(1) = B(1)  
DO  I = N-1, 1, -1
+
DO  I = 2, N-1
S = Y(I)
+
S = B(I)
DO  J = N, I+1, -1
+
DO  J = 1, I-1
S = S - DPROD(U(I,J), X(J))
+
S = S - DPROD(L(I,J), Y(J))
END DO
+
END DO
X(I) = S / U(I,I)
 
END DO
 
 
END DO
 
END DO
 
</source>
 
</source>
 
При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности.
 
При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности.
 +
 +
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
=== Результаты прогонов ===
 +
=== Выводы для классов архитектур ===
  
 
== Литература ==
 
== Литература ==
  
# В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. М.: Наука, 1984, стр. 182.
+
[[En:Forward substitution]]
 +
 
 +
[[Категория:Статьи в работе]]

Текущая версия на 09:48, 18 июля 2022


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

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

Общая структура прямой подстановки с неособенной нижней треугольной матрицей, тем не менее, практически полностью совпадает со структурой обратной подстановки. Поэтому здесь мы рассмотрим случай, когда матрица [math]L[/math], как полученная из разложения типа Гаусса, имеет единичные диагональные элементы.

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

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

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

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

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

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

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

Вычислительное ядро прямой подстановки можно составить из множественных (всего их [math]n-1[/math]) вычислений скалярных произведений строк матрицы [math]L[/math] на уже вычисленную часть вектора [math]y[/math]:

[math] \sum_{j = 1}^{i-1} l_{ij} y_{j} [/math]

в режиме накопления или без него, в зависимости от требований задачи, с их последующим вычитанием из компоненты вектора [math]b[/math] и деления на диагональный элемент матрицы [math]L[/math]. В отечественных реализациях, даже в последовательных, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм типа

[math] b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} [/math]

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

[math] b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} [/math]

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

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

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

[math] b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j} [/math]

в режиме накопления или без него.

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

Чтобы понять последовательность исполнения, перепишем формулы метода так:

1. [math]y_{1} = b_{1}[/math]

Далее для всех [math]i[/math] от [math]2[/math] до [math]n[/math] по возрастанию выполняются

2. [math]y_{i} = b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j}[/math]

Особо отметим, что вычисления сумм вида [math]b_{i} - \sum_{j = 1}^{i-1} l_{ij} y_{j}[/math] производят в режиме накопления вычитанием из [math]b_{i}[/math] произведений [math]l_{ij} y_{j}[/math] для [math]j[/math] от [math]1[/math] до [math]i-1[/math], c возрастанием [math]j[/math]. Другие порядки выполнения суммирования приводят к резкому ухудшению параллельных свойств алгоритма.

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

Для прямой подстановки порядка n в последовательном (наиболее быстром) варианте требуется:

  • по [math]\frac{n^2-n}{2}[/math] умножений и сложений (вычитаний).

При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает затраты во времени, требуемом для выполнения прямой подстановки.

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

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

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

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

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

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

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

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

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

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

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

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

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

Для прямой подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 Результаты прогонов

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

3 Литература

  1. В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. М.: Наука, 1984, стр. 182.