Классическая монотонная прогонка, повторный вариант

Повторная прогонка для трёхдиагональной матрицы, точечный вариант
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ 5n-4 }[/math]
Объём входных данных	[math]\displaystyle{ 4n-2 }[/math]
Объём выходных данных	[math]\displaystyle{ n }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ 5n-4 }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ 1 }[/math]

Основные авторы описания: А.В.Фролов.

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

Рисунок 1. Граф алгоритма повторной прогонки при n=7 без отображения входных и выходных данных. / - деление, f - операция a+bc.

Общее описание алгоритма

Повторная прогонка - один из вариантов метода исключения неизвестных, применяемого к решению СЛАУ вида [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]

Часто, однако, при изложении сути методов решения трёхдиагональных СЛАУ элементы правой части и матрицы системы обозначают и нумеруют по-другому, например СЛАУ может иметь вид ([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].

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

Рисунок 2. Граф алгоритма повторной прогонки с заранее выполненным предвычислением обратных чисел при n=7 без отображения входных и выходных данных. m - умножение на предвычисленное обратное число, f - операция a+bc.

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

В приведённых выше обозначениях в повторной прогонке сначала выполняют её прямой ход - вычисляют коэффициенты

[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].

В литературе указывается, что данные формулы эквивалентны использованию предвычисленного при полной прогонке одного из вариантов [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

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

Как видно из анализа графа алгоритма, его (без существенных модификаций) практически невозможно распараллелить. Поэтому есть два способа использования прогонок для параллельных вычислительных систем: либо разбивать задачу, где используются прогонки, так, чтобы их было достаточно много, например, так, чтобы на каждую из прогонок приходился 1 процессор (1 ядро), либо использовать вместо прогонки её параллельные варианты (циклическую редукцию, последовательно-параллельные варианты и т.п.).

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

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

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