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

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

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

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

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

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

Повторная прогонка - один из вариантов метода исключения неизвестных, применяемого к решению СЛАУ вида $Ax = b$, где как минимум один раз уже была решена прогонкой другая СЛАУ с той же матрицей

$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}$

Часто, однако, при изложении сути методов решения трёхдиагональных СЛАУ элементы правой части и матрицы системы обозначают и нумеруют по-другому, например СЛАУ может иметь вид ($N=n-1$)

$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}$

или, если записывать отдельно по уравнениям, то

$c_{0} y_{0} - b_{0} y_{1} = f_{0}$,

$-a_{i} y_{i-1} + c_{i} y_{i} - b_{i} y_{i+1} = f_{i}, 1 \le i \le N-1$,

$-a_{N} y_{N-1} + c_{N} y_{N} = f_{N}$.

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

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

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

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

$\beta_{1} = f_{0}/c_{0}$,

$\beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i})$, $\quad i = 1, 2, \cdots , N$.

(при этом отношения $1/c_{0}$ и $1/(c_{i}-a_{i}\alpha_{i})$ уже предвычислены при решении первой системы первой прогонкой) после чего вычисляют решение с помощью обратного хода

$y_{N} = \beta_{N+1}$,

$y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}$, $\quad i = N-1, N-2, \cdots , 1, 0$.

В литературе указывается, что данные формулы эквивалентны использованию предвычисленного при полной прогонке одного из вариантов $LU$-разложения матрицы системы для последующего решения двухдиагональных систем посредством прямой и обратной подстановок.

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

Вычислительное ядро алгоритма можно считать состоящим из двух частей - прямого и обратного хода прогонки. Как в прямом ходе, так и в обратном, вычислительное ядро составляют последовательности операций умножения (в прямом ходе их две за шаг; в принципе, можно воспользоваться ассоциативностью и вычислять сразу оба произведения, но "промежуточный" результат нужен для обратного хода) и сложения/вычитания.

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

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

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

Последовательность исполнения метода следующая:

1. Инициализируется прямой ход прогонки:

$\beta_{1} = f_{0}(1/c_{0})$

2. Последовательно для всех $i$ от $1$ до $N-1$ выполняются формулы прямого хода:

$\beta_{i+1} = (f_{i}+a_{i}\beta_{i})(1/(c_{i}-a_{i}\alpha_{i}))$.

3. Инициализируется обратный ход прогонки:

$y_{N} = (f_{N}+a_{N}\beta_{N})(1/(c_{N}-a_{N}\alpha_{N}))$

4. Последовательно для всех $i$ с убыванием от $N-1$ до $0$ выполняются формулы обратного хода: $y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}$.

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

Для выполнения повторной прогонки в трёхдиагональной СЛАУ из $n$ уравнений с $n$ неизвестными в последовательном варианте требуется:

• $2n-2$ сложений/вычитаний,
• $3n-2$ умножений.

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

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

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

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

Для выполнения повторной прогонки в трёхдиагональной СЛАУ из $n$ уравнений с $n$ неизвестными в параллельном варианте требуется последовательно выполнить следующие ярусы:

• $3n - 2$ ярусов умножений и $2n - 2$ сложений/вычитаний (в ярусах - по $1$ операции).

Таким образом, при классификации по высоте ЯПФ, прогонка относится к алгоритмам со сложностью $O(n)$. При классификации по ширине ЯПФ её сложность будет равна $1$ (чисто последовательный алгоритм).

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

Входные данные: Предварительно обработанная полной прогонкой трёхдиагональная матрица $A$ (элементы $a_{ij}$), вектор $b$ (элементы $b_{i}$).

Объём входных данных: $4n-2$.

Выходные данные: вектор $x$ (элементы $x_{i}$).

Объём выходных данных: $n$.

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

Соотношение последовательной и параллельной сложности, как хорошо видно, равно $1$.

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

Алгоритм в рамках выбранной версии полностью детерминирован.

Обычно прогонка используется для решения СЛАУ с диагональным преобладанием. В этом случае гарантируется устойчивость алгоритма.

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

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

В зависимости от способа хранения матрицы СЛАУ (в виде одного массива с $3$ строками или в виде $3$ разных массивов) и способа хранения вычисляемых коэффициентов (на месте уже использованных элементов матрицы либо отдельно) возможны различные реализации алгоритма.

Приведем пример подпрограммы на языке 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.4 Выводы для классов архитектур

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