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

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

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

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

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

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

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

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

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

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

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

• [math]3n - 2[/math] ярусов умножений и [math]2n - 2[/math] сложений/вычитаний (в [math]n-1[/math] ярусах - по [math]2[/math] операции, в [math]n-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 Локальность данных и вычислений

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

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

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

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

2.4 Масштабируемость алгоритма и его реализации

О масштабируемости самой прогонки, как полностью непараллельного алгоритма, говорить не имеет смысла. Однако необходимо отметить, что анализ масштабируемости параллельных вариантов прогонки должен проводиться относительно однопроцессорной реализации описанного классического варианта прогонки, а не относительно однопроцессорных расчетов для её параллельных вариантов.

2.5 Динамические характеристики и эффективность реализации алгоритма

В силу существенно последовательной природы алгоритма и его избыточной локальности, исследование его динамических характеристик представляется малоценным.

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

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

2.7 Существующие реализации алгоритма

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