Классическая монотонная прогонка, повторный вариант
| Повторная прогонка для трёхдиагональной матрицы, точечный вариант | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ 5n-4 }[/math] |
| Объём входных данных | [math]\displaystyle{ 4n-2 }[/math] |
| Объём выходных данных | [math]\displaystyle{ n }[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]\displaystyle{ 5n-4 }[/math] |
| Ширина ярусно-параллельной формы | [math]\displaystyle{ 1 }[/math] |
Основные авторы описания: А.В.Фролов
Свойства и структура алгоритма
Общее описание алгоритма
Повторная прогонка - один из вариантов метода исключения неизвестных, применяемого к решению СЛАУ[1][2] вида [math]\displaystyle{ Ax = b }[/math], где как минимум один раз уже была решена прогонкой другая СЛАУ с той же матрицей
- [math]\displaystyle{ A = \begin{bmatrix} a_{11} & a_{12} & 0 & \cdots & \cdots & 0 \\ a_{21} & a_{22} & a_{23}& \cdots & \cdots & 0 \\ 0 & a_{32} & a_{33} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & a_{n-1 n-2} & a_{n-1 n-1} & a_{n-1 n} \\ 0 & \cdots & \cdots & 0 & a_{n n-1} & a_{n n} \\ \end{bmatrix}, x = \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \\ \end{bmatrix}, b = \begin{bmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \\ \end{bmatrix} }[/math]
Часто, однако, при изложении сути методов решения трёхдиагональных СЛАУ[3] элементы правой части и матрицы системы обозначают и нумеруют по-другому, например СЛАУ может иметь вид ([math]\displaystyle{ N=n-1 }[/math])
- [math]\displaystyle{ A = \begin{bmatrix} c_{0} & -b_{0} & 0 & \cdots & \cdots & 0 \\ -a_{1} & c_{1} & -b_{1} & \cdots & \cdots & 0 \\ 0 & -a_{2} & c_{2} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & -a_{N-1} & c_{N-1} & -b_{N-1} \\ 0 & \cdots & \cdots & 0 & -a_{N} & c_{N} \\ \end{bmatrix}\begin{bmatrix} y_{0} \\ y_{1} \\ \vdots \\ y_{N} \\ \end{bmatrix} = \begin{bmatrix} f_{0} \\ f_{1} \\ \vdots \\ f_{N} \\ \end{bmatrix} }[/math]
или, если записывать отдельно по уравнениям, то
[math]\displaystyle{ c_{0} y_{0} - b_{0} y_{1} = f_{0} }[/math],
[math]\displaystyle{ -a_{i} y_{i-1} + c_{i} y_{i} - b_{i} y_{i+1} = f_{i}, 1 \le i \le N-1 }[/math],
[math]\displaystyle{ -a_{N} y_{N-1} + c_{N} y_{N} = f_{N} }[/math].
Здесь рассматривается тот вариант прогонки, когда обрабатывается вся СЛАУ сверху вниз и обратно - так называемая правая прогонка. Суть метода состоит в исключении из уравнений неизвестных, сначала - сверху вниз - под диагональю, а затем - снизу вверх - над диагональю. Вариант, когда СЛАУ "проходится" наоборот, снизу вверх и обратно вниз - левая прогонка - принципиально ничем не отличается от рассматриваемого, поэтому нет смысла включать его отдельное описание.
Математическое описание алгоритма
В приведённых выше обозначениях в повторной прогонке сначала выполняют её прямой ход - вычисляют коэффициенты
[math]\displaystyle{ \beta_{1} = f_{0}/c_{0} }[/math],
[math]\displaystyle{ \beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i}) }[/math], [math]\displaystyle{ \quad i = 1, 2, \cdots , N }[/math].
(при этом отношения [math]\displaystyle{ 1/c_{0} }[/math] и [math]\displaystyle{ 1/(c_{i}-a_{i}\alpha_{i}) }[/math] уже предвычислены при решении первой системы первой прогонкой) после чего вычисляют решение с помощью обратного хода
[math]\displaystyle{ y_{N} = \beta_{N+1} }[/math],
[math]\displaystyle{ y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1} }[/math], [math]\displaystyle{ \quad i = N-1, N-2, \cdots , 1, 0 }[/math].
В литературе[3] указывается, что данные формулы эквивалентны использованию предвычисленного при полной прогонке одного из вариантов [math]\displaystyle{ LU }[/math]-разложения матрицы системы для последующего решения двухдиагональных систем посредством прямой и обратной подстановок.
Вычислительное ядро алгоритма
Вычислительное ядро алгоритма можно считать состоящим из двух частей - прямого и обратного хода прогонки. Как в прямом ходе, так и в обратном, вычислительное ядро составляют последовательности операций умножения (в прямом ходе их две за шаг; в принципе, можно воспользоваться ассоциативностью и вычислять сразу оба произведения, но "промежуточный" результат нужен для обратного хода) и сложения/вычитания.
Макроструктура алгоритма
В дополнение к возможности представления макроструктуры алгоритма как совокупности прямого и обратного хода, прямой ход также может быть разложен на две макроединицы - треугольного разложения матрицы и прямого хода решения двухдиагональной СЛАУ, которые выполняются "одновременно", т.е. параллельно друг другу. При этом решение двухдиагональной СЛАУ использует результаты треугольного разложения.
Схема реализации последовательного алгоритма
Последовательность исполнения метода следующая:
1. Инициализируется прямой ход прогонки:
[math]\displaystyle{ \beta_{1} = f_{0}(1/c_{0}) }[/math]
2. Последовательно для всех [math]\displaystyle{ i }[/math] от [math]\displaystyle{ 1 }[/math] до [math]\displaystyle{ N-1 }[/math] выполняются формулы прямого хода:
[math]\displaystyle{ \beta_{i+1} = (f_{i}+a_{i}\beta_{i})(1/(c_{i}-a_{i}\alpha_{i})) }[/math].
3. Инициализируется обратный ход прогонки:
[math]\displaystyle{ y_{N} = (f_{N}+a_{N}\beta_{N})(1/(c_{N}-a_{N}\alpha_{N})) }[/math]
4. Последовательно для всех [math]\displaystyle{ i }[/math] с убыванием от [math]\displaystyle{ N-1 }[/math] до [math]\displaystyle{ 0 }[/math] выполняются формулы обратного хода: [math]\displaystyle{ y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1} }[/math].
Последовательная сложность алгоритма
Для выполнения повторной прогонки в трёхдиагональной СЛАУ из [math]\displaystyle{ n }[/math] уравнений с [math]\displaystyle{ n }[/math] неизвестными в последовательном варианте требуется:
- [math]\displaystyle{ 2n-2 }[/math] сложений/вычитаний,
- [math]\displaystyle{ 3n-2 }[/math] умножений.
Таким образом, при классификации по последовательной сложности, повторная прогонка относится к алгоритмам с линейной сложностью.
Информационный граф
Информационный граф повторной прогонки в случае первой прогонки без предвычислений обратных чисел на рис.1. Информационный граф повторной прогонки в случае первой прогонки с предвычислениями обратных чисел - на рис.2. Как следует из анализа графа, он является полностью последовательным. Из рисунка видно, что не только математическая суть обработки элементов векторов, но даже структура графа алгоритма и направление потоков данных в нём вполне соответствуют названиям "прямой и обратный ход".
Описание ресурса параллелизма алгоритма
Для выполнения повторной прогонки в трёхдиагональной СЛАУ из [math]\displaystyle{ n }[/math] уравнений с [math]\displaystyle{ n }[/math] неизвестными в параллельном варианте требуется последовательно выполнить следующие ярусы:
- [math]\displaystyle{ 3n - 2 }[/math] ярусов умножений и [math]\displaystyle{ 2n - 2 }[/math] сложений/вычитаний (в ярусах - по [math]\displaystyle{ 1 }[/math] операции).
Таким образом, при классификации по высоте ЯПФ, прогонка относится к алгоритмам со сложностью [math]\displaystyle{ O(n) }[/math]. При классификации по ширине ЯПФ её сложность будет равна [math]\displaystyle{ 1 }[/math] (чисто последовательный алгоритм).
Входные и выходные данные алгоритма
Входные данные: Предварительно обработанная полной прогонкой трёхдиагональная матрица [math]\displaystyle{ A }[/math] (элементы [math]\displaystyle{ a_{ij} }[/math]), вектор [math]\displaystyle{ b }[/math] (элементы [math]\displaystyle{ b_{i} }[/math]).
Объём входных данных: [math]\displaystyle{ 4n-2 }[/math].
Выходные данные: вектор [math]\displaystyle{ x }[/math] (элементы [math]\displaystyle{ x_{i} }[/math]).
Объём выходных данных: [math]\displaystyle{ n }[/math].
Свойства алгоритма
Соотношение последовательной и параллельной сложности, как хорошо видно, равно [math]\displaystyle{ 1 }[/math].
При этом вычислительная мощность алгоритма как отношение числа операций к суммарному объему входных и выходных данных также является константой.
Алгоритм в рамках выбранной версии полностью детерминирован.
Обычно прогонка используется для решения СЛАУ с диагональным преобладанием. В этом случае гарантируется устойчивость алгоритма.
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
В зависимости от способа хранения матрицы СЛАУ (в виде одного массива с [math]\displaystyle{ 3 }[/math] строками или в виде [math]\displaystyle{ 3 }[/math] разных массивов) и способа хранения вычисляемых коэффициентов (на месте уже использованных элементов матрицы либо отдельно) возможны различные реализации алгоритма.
Приведем пример подпрограммы на языке Fortran, реализующей прогонку, где все элементы матрицы хранятся в одном массиве, причём соседние элементы матричной строки размещаются рядом, а вычисляемые в процессе первой прогонки коэффициенты --- на месте элементов исходной матрицы.
subroutine progmr (a,x,N)
real a(3,0:N), x(0:N)
x(0)=x(0)*a(2,0) ! beta 1
do 10 i=1,N-1
x(i)=(x(i)-a(1,i)*x(i-1))*a(2,i) ! beta i+1
10 continue
x(N)=(x(N)-a(1,N)*x(N-1))*a(2,N) ! y N
do 20 i=N-1,0,-1
x(i)=a(3,i)*x(i+1)+x(i) ! y i
20 continue
return
end
Локальность данных и вычислений
Как видно по графу алгоритма, локальность данных по пространству хорошая - все аргументы, которые необходимы для операций, вычисляются "рядом". Однако по времени локальность вычислений не столь хороша. Если данные задачи не помещаются в кэш, то вычисления в "верхнем левом углу" СЛАУ будут выполняться с постоянными промахами кэша. Отсюда может следовать одна из рекомендаций прикладникам, использующим прогонку, - нужно организовать все вычисления так, чтобы прогонки были "достаточно коротки" для помещения данных в кэш.
При этом, однако, повторная прогонка относится к таким последовательным алгоритмам, в которых локальность вычислений настолько велика, что является даже излишней[4]. Из-за того, что данные, необходимые для выполнения основных операций прогонки, вычисляются в непосредственно предшествующим им операциях, возможность использования суперскалярности вычислительных ядер процессоров практически сводится на нет, что ухудшает эффективность выполнения прогонки даже на современных однопроцессорных и одноядерных системах. Последнее показано сравнением времени исполнения монотонной повторной прогонки и встречной монотонной прогонки[5], которая благодаря наличию двух ветвей вычислений даёт выигрыш по времени около 20% в сравнении с первой даже на персональных компьютерах (см. Рисунок 2).
Возможные способы и особенности параллельной реализации алгоритма
Как видно из анализа графа алгоритма, его (без существенных модификаций) практически невозможно распараллелить. Поэтому есть два способа использования прогонок для параллельных вычислительных систем: либо разбивать задачу, где используются прогонки, так, чтобы их было достаточно много, например, так, чтобы на каждую из прогонок приходился 1 процессор (1 ядро), либо использовать вместо прогонки её параллельные варианты (циклическую редукцию, последовательно-параллельные варианты и т.п.).
Масштабируемость алгоритма и его реализации
О масштабируемости самой прогонки, как полностью непараллельного алгоритма, говорить не имеет смысла. Однако необходимо отметить, что анализ масштабируемости параллельных вариантов прогонки должен проводиться относительно однопроцессорной реализации описанного классического варианта прогонки, а не относительно однопроцессорных расчетов для её параллельных вариантов.
Динамические характеристики и эффективность реализации алгоритма
В силу существенно последовательной природы алгоритма и его избыточной локальности, исследование его динамических характеристик представляется малоценным.
Выводы для классов архитектур
Повторная монотонная прогонка - метод для архитектуры классического, фон-неймановского типа, устаревший даже для эффективной загрузки одноядерных систем, поддерживающих суперскалярность, где проигрывает повторной встречной прогонке. Для распараллеливания решения СЛАУ с трёхдиагональной матрицей следует использовать какой-либо её параллельный заменитель, например, наиболее распространённую циклическую редукцию или уступающий ей по критическому пути графа, но имеющий более регулярную структуру графа новый последовательно-параллельный вариант.
Существующие реализации алгоритма
Алгоритм прогонки является настолько простым, что, несмотря на его "дежурное" присутствие в стандартных пакетах программ решения задач линейной алгебры, большинство использующих его исследователей-прикладников часто пишут соответствующий фрагмент программы самостоятельно. Однако надо отметить, что, как правило, метод прогонки в пакетах реализуют не в его "чистом виде", а в виде пары "разложение на две двухдиагональные матрицы" и "решение СЛАУ с предварительно факторизованной трёхдиагональной матрицей". Так обстоит дело, например, в пакете SCALAPACK и его предшественниках.
Литература
- ↑ Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
- ↑ Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.
- ↑ 3,0 3,1 Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978.
- ↑ Фролов А.В., Антонов А.С., Воеводин Вл.В., Теплов А.М. Сопоставление разных методов решения одной задачи по методике проекта Algowiki // Параллельные вычислительные технологии (ПаВТ'2016): труды международной научной конференции (г. Архангельск, 28 марта - 1 апреля 2016 г.). Челябинск: Издательский центр ЮУрГУ, 2016. С. 347-360.
- ↑ Фролов Н.А., Фролов А.В. Экспериментальные исследования влияния степени локальности алгоритмов на их быстродействие на примере решения трёхдиагональных СЛАУ // Труды 59й научной конференции МФТИ (21–26 ноября 2016 г., гг. Москва-Долгопрудный).