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

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

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

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

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

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

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

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

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. Информационный граф повторной прогонки в случае первой прогонки с предвычислениями обратных чисел - на рис.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 Локальность данных и вычислений

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

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

Рисунок 2. Отношение времени исполнения повторной встречной и монотонных прогонок. По оси асбсцисс отложено [math]m[/math] (порядок СЛАУ равен [math]2^m - 1[/math]), по оси ординат - отношение времени, затраченного на исполнение встречной повторной прогонки, к времени, затраченному на выполнение монотонной повторной прогонки. Цветами выделены графики для разных систем (Pentium D + Windows XP, Core i5 + Windows 8, Core i7 + Windows 7), на всех их использована одна и та же сборка кода компилятором GNU Fortran.

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

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

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

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

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

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

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

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

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

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