Классическая монотонная прогонка, повторный вариант
| Повторная прогонка для трёхдиагональной матрицы, точечный вариант | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]5n-4[/math] |
| Объём входных данных | [math]4n-2[/math] |
| Объём выходных данных | [math]n[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]5n-4[/math] |
| Ширина ярусно-параллельной формы | [math]1[/math] |
Основные авторы описания: А.В.Фролов.
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Описание ресурса параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Повторная прогонка - один из вариантов метода исключения неизвестных, применяемого к решению СЛАУ[1][2] вида [math]Ax = b[/math], где как минимум один раз уже была решена прогонкой другая СЛАУ с той же матрицей
- [math] 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]N=n-1[/math])
- [math] 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]c_{0} y_{0} - b_{0} y_{1} = f_{0}[/math],
[math]-a_{i} y_{i-1} + c_{i} y_{i} - b_{i} y_{i+1} = f_{i}, 1 \le i \le N-1[/math],
[math]-a_{N} y_{N-1} + c_{N} y_{N} = f_{N}[/math].
Здесь рассматривается тот вариант прогонки, когда обрабатывается вся СЛАУ сверху вниз и обратно - так называемая правая прогонка. Суть метода состоит в исключении из уравнений неизвестных, сначала - сверху вниз - под диагональю, а затем - снизу вверх - над диагональю. Вариант, когда СЛАУ "проходится" наоборот, снизу вверх и обратно вниз - левая прогонка - принципиально ничем не отличается от рассматриваемого, поэтому нет смысла включать его отдельное описание.
1.2 Математическое описание алгоритма
В приведённых выше обозначениях в повторной прогонке сначала выполняют её прямой ход - вычисляют коэффициенты
[math]\beta_{1} = f_{0}/c_{0}[/math],
[math]\beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i})[/math], [math]\quad i = 1, 2, \cdots , N[/math].
(при этом отношения [math]1/c_{0}[/math] и [math]1/(c_{i}-a_{i}\alpha_{i})[/math] уже предвычислены при решении первой системы первой прогонкой) после чего вычисляют решение с помощью обратного хода
[math]y_{N} = \beta_{N+1}[/math],
[math]y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}[/math], [math]\quad i = N-1, N-2, \cdots , 1, 0[/math].
В литературе[3] указывается, что данные формулы эквивалентны использованию предвычисленного при полной прогонке одного из вариантов [math]LU[/math]-разложения матрицы системы для последующего решения двухдиагональных систем посредством прямой и обратной подстановок.
1.3 Вычислительное ядро алгоритма
Вычислительное ядро алгоритма можно считать состоящим из двух частей - прямого и обратного хода прогонки. Как в прямом ходе, так и в обратном, вычислительное ядро составляют последовательности операций умножения (в прямом ходе их две за шаг; в принципе, можно воспользоваться ассоциативностью и вычислять сразу оба произведения, но "промежуточный" результат нужен для обратного хода) и сложения/вычитания.
1.4 Макроструктура алгоритма
В дополнение к возможности представления макроструктуры алгоритма как совокупности прямого и обратного хода, прямой ход также может быть разложен на две макроединицы - треугольного разложения матрицы и прямого хода решения двухдиагональной СЛАУ, которые выполняются "одновременно", т.е. параллельно друг другу. При этом решение двухдиагональной СЛАУ использует результаты треугольного разложения.
1.5 Схема реализации последовательного алгоритма
Последовательность исполнения метода следующая:
1. Инициализируется прямой ход прогонки:
[math]\beta_{1} = f_{0}(1/c_{0}) [/math]
2. Последовательно для всех [math]i[/math] от [math]1[/math] до [math]N-1[/math] выполняются формулы прямого хода:
[math]\beta_{i+1} = (f_{i}+a_{i}\beta_{i})(1/(c_{i}-a_{i}\alpha_{i}))[/math].
3. Инициализируется обратный ход прогонки:
[math]y_{N} = (f_{N}+a_{N}\beta_{N})(1/(c_{N}-a_{N}\alpha_{N}))[/math]
4. Последовательно для всех [math]i[/math] с убыванием от [math]N-1[/math] до [math]0[/math] выполняются формулы обратного хода: [math]y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}[/math].
1.6 Последовательная сложность алгоритма
Для выполнения повторной прогонки в трёхдиагональной СЛАУ из [math]n[/math] уравнений с [math]n[/math] неизвестными в последовательном варианте требуется:
- [math]2n-2[/math] сложений/вычитаний,
- [math]3n-2[/math] умножений.
Таким образом, при классификации по последовательной сложности, повторная прогонка относится к алгоритмам с линейной сложностью.
1.7 Информационный граф
Информационный граф повторной прогонки в случае первой прогонки без предвычислений обратных чисел на рис.1. Информационный граф повторной прогонки в случае первой прогонки с предвычислениями обратных чисел - на рис.2. Как следует из анализа графа, он является полностью последовательным. Из рисунка видно, что не только математическая суть обработки элементов векторов, но даже структура графа алгоритма и направление потоков данных в нём вполне соответствуют названиям "прямой и обратный ход".
1.8 Описание ресурса параллелизма алгоритма
Для выполнения повторной прогонки в трёхдиагональной СЛАУ из [math]n[/math] уравнений с [math]n[/math] неизвестными в параллельном варианте требуется последовательно выполнить следующие ярусы:
- [math]3n - 2[/math] ярусов умножений и [math]2n - 2[/math] сложений/вычитаний (в ярусах - по [math]1[/math] операции).
Таким образом, при классификации по высоте ЯПФ, прогонка относится к алгоритмам со сложностью [math]O(n)[/math]. При классификации по ширине ЯПФ её сложность будет равна [math]1[/math] (чисто последовательный алгоритм).
1.9 Входные и выходные данные алгоритма
Входные данные: Предварительно обработанная полной прогонкой трёхдиагональная матрица [math]A[/math] (элементы [math]a_{ij}[/math]), вектор [math]b[/math] (элементы [math]b_{i}[/math]).
Объём входных данных: [math]4n-2[/math].
Выходные данные: вектор [math]x[/math] (элементы [math]x_{i}[/math]).
Объём выходных данных: [math]n[/math].
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности, как хорошо видно, равно [math]1[/math].
При этом вычислительная мощность алгоритма как отношение числа операций к суммарному объему входных и выходных данных также является константой.
Алгоритм в рамках выбранной версии полностью детерминирован.
Обычно прогонка используется для решения СЛАУ с диагональным преобладанием. В этом случае гарантируется устойчивость алгоритма.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
В зависимости от способа хранения матрицы СЛАУ (в виде одного массива с [math]3[/math] строками или в виде [math]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
2.2 Возможные способы и особенности параллельной реализации алгоритма
Как видно из анализа графа алгоритма, его (без существенных модификаций) практически невозможно распараллелить. Поэтому есть два способа использования прогонок для параллельных вычислительных систем: либо разбивать задачу, где используются прогонки, так, чтобы их было достаточно много, например, так, чтобы на каждую из прогонок приходился 1 процессор (1 ядро), либо использовать вместо прогонки её параллельные варианты (циклическую редукцию, последовательно-параллельные варианты и т.п.).
Алгоритм прогонки является настолько простым, что, несмотря на его "дежурное" присутствие в стандартных пакетах программ решения задач линейной алгебры, большинство использующих его исследователей-прикладников часто пишут соответствующий фрагмент программы самостоятельно. Однако надо отметить, что, как правило, метод прогонки в пакетах реализуют не в его "чистом виде", а в виде пары "разложение на две двухдиагональные матрицы" и "решение СЛАУ с предварительно факторизованной трёхдиагональной матрицей". Так обстоит дело, например, в пакете SCALAPACK и его предшественниках.
2.3 Результаты прогонов
2.4 Выводы для классов архитектур
Повторная монотонная прогонка - метод для архитектуры классического, фон-неймановского типа, устаревший даже для эффективной загрузки одноядерных систем, поддерживающих суперскалярность, где проигрывает повторной встречной прогонке. Для распараллеливания решения СЛАУ с трёхдиагональной матрицей следует использовать какой-либо её параллельный заменитель, например, наиболее распространённую циклическую редукцию или уступающий ей по критическому пути графа, но имеющий более регулярную структуру графа новый последовательно-параллельный вариант.
