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

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

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][досмотренная версия]
 
(не показано 20 промежуточных версий 6 участников)
Строка 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>).
Строка 49: Строка 60:
 
в режиме накопления или без него.
 
в режиме накопления или без него.
  
=== Описание схемы реализации последовательного алгоритма ===
+
=== Схема реализации последовательного алгоритма ===
  
 
Чтобы понять последовательность исполнения, перепишем формулы метода так:
 
Чтобы понять последовательность исполнения, перепишем формулы метода так:
Строка 75: Строка 86:
 
Опишем [[глоссарий#Граф алгоритма|граф алгоритма]] как аналитически, так и в виде рисунка.  
 
Опишем [[глоссарий#Граф алгоритма|граф алгоритма]] как аналитически, так и в виде рисунка.  
  
[[Файл:DirectL.png|450px|thumb|left|Прямая подстановка]]
+
[[Файл:DirectL.png|450px|thumb|left|Рисунок 1. Прямая подстановка]]
  
 
Граф алгоритма прямой подстановки состоит из одной группы вершин, расположенной в целочисленных узлах двумерной области, соответствующая ей операция  <math>a-bc</math>.  
 
Граф алгоритма прямой подстановки состоит из одной группы вершин, расположенной в целочисленных узлах двумерной области, соответствующая ей операция  <math>a-bc</math>.  
Строка 93: Строка 104:
  
 
Результат срабатывания операции является:
 
Результат срабатывания операции является:
** при <math>j < i-1</math> - ''промежуточным данным'' алгоритма;
+
* при <math>j < i-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>. Здесь вершины обозначены зелёным цветом. Подача входных данных из вектора <math>b</math>, идущая в вершины левого столбца, и матрицы <math>L</math>, идущая во все вершины, на рисунке не представлена.
  
=== Описание ресурса параллелизма алгоритма ===
+
=== Ресурс параллелизма алгоритма ===
  
 
Для прямой подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:
 
Для прямой подстановки порядка n в параллельном варианте требуется последовательно выполнить следующие ярусы:
Строка 108: Строка 119:
 
При классификации по высоте ЯПФ, таким образом, метод прямой подстановки относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''линейной''.
 
При классификации по высоте ЯПФ, таким образом, метод прямой подстановки относится к алгоритмам ''с линейной сложностью''. При классификации по ширине ЯПФ его сложность также будет ''линейной''.
  
=== Описание входных и выходных данных ===
+
=== Входные и выходные данные алгоритма ===
  
'''Входные данные''': левая треугольная матрица <math>L</math> (элементы <math>l_{ij}</math>), вектор правой части <math>b</math> (элементы <math>b_{i}</math>).
+
'''Входные данные''': нижняя треугольная матрица <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>\frac{n (n + 1)}{2}</math> (в силу треугольности и единичности диагональных элементов достаточно хранить только поддиагональные элементы матрицы <math>L</math>).  
Строка 116: Строка 127:
 
'''Выходные данные''': вектор решения <math>y</math> (элементы <math>y_{i}</math>).
 
'''Выходные данные''': вектор решения <math>y</math> (элементы <math>y_{i}</math>).
  
'''Объём выходных данных''': <math>n</math>.
+
'''Объём выходных данных''': <math>n~.</math>
  
 
=== Свойства алгоритма ===
 
=== Свойства алгоритма ===
Строка 126: Строка 137:
 
При этом алгоритм прямой подстановки полностью детерминирован. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и делает параллельную сложность квадратической.
 
При этом алгоритм прямой подстановки полностью детерминирован. Использование другого порядка выполнения ассоциативных операций в данной версии нами не рассматривается, поскольку в корне меняет структуру алгоритма и делает параллельную сложность квадратической.
  
== Программная реализация ==
+
== Программная реализация алгоритма ==
  
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
Строка 143: Строка 154:
 
При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности.
 
При этом для реализации режима накопления переменная <math>S</math> должна быть двойной точности.
  
=== Описание локальности данных и вычислений ===
+
=== Возможные способы и особенности параллельной реализации алгоритма ===
==== Описание локальности алгоритма ====
+
=== Результаты прогонов ===
==== Описание локальности реализации алгоритма ====
 
===== Описание структуры обращений в память и качественная оценка локальности =====
 
 
 
[[file:gauss forward 1.PNG|thumb|center|700px|Рисунок 12.1. Прямой ход метода Гаусса решения СЛАУ. Общий профиль обращений в память]]
 
 
 
На рисунке 12.1 представлен профиль обращений в память для прямого хода метода Гаусса решения СЛАУ. Из исходного кода, как и из рисунка, видно, что в программе задействовано 2 основных массива: профиль обращений к элементам первого выделен зеленым (фрагмент 1); все остальные обращения происходят к элементам второго массива. Видно, что общий профиль всей программы устроен достаточно сложно. Основное, что можно заметить из данного рисунка, – это итерационная природа профиля, причем каждая следующая итерация содержит меньше обращений. Это особенно хорошо видно по фрагменту  1.
 
 
 
Далее, как известно, в прямом ходе Гаусса на каждой итерации отбрасывается из рассмотрения одна из строк СЛАУ. Это отчетливо видно и на рис. 12.1: на каждой итерации (одна из них выделена как фрагмент 2) выполняется серия обращений, связанных с отбрасываемой строкой (фрагмент 3), после чего обращения к этим элементам больше не выполняются. Также можно заметить, что фрагменты, подобные фрагменту 3, на каждой последующей итерации становятся все меньше, как по числу обращений, так и по числу элементов, к которым происходят обращения.
 
 
 
Однако на основе подобной общей картины очень сложно дать даже качественные оценки локальности. Перейдем к более подробному рассмотрению фрагмента 1, соответствующего обращениям к одному из двух массивов.
 
 
 
[[file:gauss forward 2.PNG|thumb|center|500px|Рисунок 12.2. Фрагмент 1 (профиль обращений к первому массиву)]]
 
 
 
На рис. 12.2 показан отдельно фрагмент 1, выделенный на рис. 12.1 зеленым. Наблюдается довольно интересная картина: исходя из рис. 12.1, складывалось впечатление, что данный профиль образуется чередованием итераций двух разных циклов, причем итерация первого цикла заметно короче итераций второго цикла. На самом же деле, из рис. 12.2 можно увидеть, что профиль состоит из однотипных итераций, в рамках которых производится последовательный перебор элементов массива. При этом на каждой следующей итерации отбрасывается элемент с минимальным индексом. Иллюзия наличия двух разных циклов обусловлена исключительно числом обращений к другому массиву – когда число таких обращений велико, итерация данного фрагмента кажется длиннее. При этом во фрагменте 1 (как видно на рис. 12.2) число обращений всего лишь около 1200, тогда как весь профиль обращений состоит из 34 тысяч. То есть на долю обращений к первому массиву приходится всего 3.5% всех обращений в программе.
 
 
 
Поскольку основа данного профиля – последовательный перебор элементов массива, можно сказать, что данный фрагмент обладает высокой пространственной локальностью. А поскольку итераций с перебором элементов достаточно много, также видно, что и временна́я локальность достаточно высока.
 
 
 
Теперь перейдем к рассмотрению профиля обращений ко второму массиву, который значительно больше и сложнее первого профиля. На рис. 12.3 отдельно представлен фрагмент 2, выделенный на рис. 12.1. Такой фрагмент устроен из четырех частей, отмеченных на рисунке оранжевым.
 
 
 
[[file:gauss forward 3.PNG|thumb|center|700px|Рисунок 12.3. Фрагмент 2 (итерация профиля обращений ко второму массиву)]]
 
 
 
Часть 1 очень похожа на повторяющийся последовательный перебор одного и того же небольшого числа элементов в цикле. При более близком рассмотрении можно убедиться, что это действительно так. На рис. 12.4 показано приближение фрагмента 2, выделенное на рис. 12.3 синим цветом. Также на этом рисунке видно, что часть 2 представляет собой практически полное подобие последовательного перебора всех элементов данного массива. Таким образом, обе эти части обладают высокой пространственной локальностью (из-за последовательного перебора). При этом часть 1 обладает также высокой временной локальностью, поскольку перебор повторяется в цикле, в отличие от части 2, где к каждому элементу выполняется только одно обращение.
 
 
 
Часть 3 состоит всего из нескольких обращений и представляет собой перебор всех элементов массива, выполненных с достаточно большим шагом, гораздо большим, чем в части 2. В данном случае, из-за достаточно большого шага,  наблюдается невысокая пространственная локальность и отсутствие временно́й локальности, так как перебор выполняется единожды.
 
 
 
Часть 4 представляет собой последовательный перебор части массива, по размеру соответствующей числу отброшенных из рассмотрения элементов; чем ближе к концу профиля, тем больше обращений расположено в части 4. Эта часть обладает теми же свойствами локальности, что и часть 2. Отличие заключается в числе обращений, однако в данном случае это не оказывает сильное влияние на локальность.
 
 
 
При рассмотрении фрагмента 2 в целом основное влияние оказывают части 1 и 2, поэтому можно утверждать, что это фрагмент обладает достаточно высокой и пространственной, и временно́й локальностью. Однако на последующих итерациях (см. рис. 12.1) локальность в целом будет постепенно снижаться, поскольку на каждой итерации отбрасывается из рассмотрения некоторая часть элементов.
 
 
 
[[file:gauss forward 4.PNG|thumb|center|500px|Рисунок 12.4. Приближение фрагмента 2]]
 
 
 
===== Количественная оценка локальности =====
 
 
 
Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.
 
 
 
На рисунке 12.5 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что отмеченный синий столбец, соответствующий прямому ходу метода Гаусса, расположен в правой части и показывает достаточно высокую производительность, немного ниже производительности теста Линпак. Также отметим, что значение daps для данной задачи практически совпадает со значением daps для всего метода Гаусса (столбец «gauss»), объединяющего как прямой, так и обратный ход. Это связано с тем, что обратный ход состоит из гораздо меньшего числа обращений, нежели прямой ход, поэтому оказывает меньшее влияние.
 
 
 
[[file:gauss forward daps ru.PNG|thumb|center|700px|Рисунок 12.5. Сравнение значений оценки daps]]
 
 
 
Вторая характеристика – cvg – предназначена для получения более машинно-независимой оценки локальности. Она определяет, насколько часто в программе необходимо подтягивать данные в кэш-память. Соответственно, чем меньше значение cvg, тем реже это нужно делать, тем лучше локальность.
 
 
 
На рисунке 12.6 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). Можно увидеть, что, согласно данной оценке, как и в случае с оценкой daps, прямой ход метода Гаусса обладает высокой локальностью. Причем согласно данной оценке тест Линпак показывает практически те же результаты. Отметим, что снова значение оценки для отдельно прямого хода и совокупности прямого и обратного ходов (столбец «gauss») совпадают.
 
 
 
[[file:gauss forward cvg.PNG|thumb|center|700px|Рисунок 12.6. Сравнение значений оценки cvg]]
 
 
 
=== Возможные способы и особенности реализации параллельного алгоритма ===
 
=== Масштабируемость алгоритма и его реализации ===
 
==== Описание масштабируемости алгоритма ====
 
==== Описание масштабируемости реализации алгоритма ====
 
[[file:Масштабируемость параллельной реализации метода Гаусса Производительность.png|thumb|center|700px|Рисунок 1. Масштабируемость параллельной реализации метода Гаусса (прямой ход) Максимальная производительность. ]]
 
[[file:Масштабируемость метод Гаусса.png|thumb|center|700px|Рисунок 2. Масштабируемость параллельной реализации метода Гаусса (прямой ход) Максимальная эффективность. ]]
 
Набор изменяемых параметров запуска реализации алгоритма и границы значений параметров алгоритма:
 
* число процессоров [4 : 256]
 
* размер матрицы [1024 : 5120]
 
Эффективность выполнения реализации алгоритма
 
* Минимальная эффективность 0,11 %
 
* Максимальная эффективность 6.65%
 
Оценка масштабируемости
 
* По числу процессов: -0.000101 – при увеличении числа процессов эффективность в целом уменьшается по рассматриваемой области, хотя и менее интенсивно, чем при увеличении числа процессов. Учитывая разницу между максимальным и минимальным значением эффективности в 6,5% можно сделать вывод, что на рассмотренной области снижение эффективности, скорее всего, достаточно равномерное. Снижение эффективности объясняется тем, что при росте вычислительной сложности существенно возрастают объемы передаваемых данных. Могут присутствовать области возрастания эффективности, на всех рассмотренных размерах матрицы. Это может объясняться увеличением вычислительной сложности задачи, что при постоянных накладных расходах на организацию параллельного взаимодействия приводит к общему увеличению эффективности работы.
 
* По размеру задачи:-0.00674– при увеличении размера задачи эффективность в целом уменьшается по рассматриваемой области. Это может свидетельствовать о том, что при увеличении размера задачи возрастают и накладные расходы на организацию параллельного взаимодействия. Так как интенсивность снижения эффективности по этому направлению выше, чем по направлению числа процессов, скорее всего падение интенсивности значительно при малом числе процессов более интенсивное, чем в присутствующих областях возрастания эффективности при большом числе процессов.
 
* По двум направлениям: -6.621e-05 – при рассмотрении увеличения, как вычислительной сложности, так и числа процессов по всей рассмотренной области значений уменьшается, однако интенсивность уменьшения эффективности очень мала. В совокупности с тем фактом, что разница между максимальной и минимальной эффективностью на рассмотренной области значений параметров составляет почти 6,5 % говорит о том, что если на поверхности присутствуют области с очень интенсивным изменением эффективности, но очень малые по площади. На остальной поверхности изменения эффективности незначительны и находятся на приблизительно одном и том же уровне.
 
 
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
 
=== Выводы для классов архитектур ===
 
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
 
  
 
== Литература ==
 
== Литература ==
  
# В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. М.: Наука, 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.