Повторная прогонка, точечный вариант

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

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

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

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

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

Повторная прогонка - часть метода прогонки в приложении к решению СЛАУ вида

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

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

Рисунок 2. Граф алгоритма повторной прогонки с предвычислением обратных чисел при n=8 без отображения входных и выходных данных. 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]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}/c_{0}[/math],

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

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

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

[math] y_{N} = (f_{N}+a_{N}\beta_{N})/(c_{N}-a_{N}\alpha_{N}) [/math]

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

[math] y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}. [/math]

Для быстрейшей работы повторной прогонки все деления в формулах заменяются умножением на предвычисленные обратные числа ([math]1/c_{0}[/math] и [math]1/(c_{i}-a_{i}\alpha_{i})[/math].

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

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

• [math]2n-2[/math] сложений,
• [math]3n-2[/math] умножений.

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

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

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

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

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

• [math]n[/math] ярусов умножений (в каждом из ярусов по одному умножению),
• по [math]2n - 2[/math] ярусов умножений и сложений (все ярусы - по одной операции).

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

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

Соотношение последовательной и параллельной сложности, как хорошо видно, является константой (причём это единица).

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

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

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

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

В зависимости от нужд вычислений, возможны как разные способы хранения разложения матрицы СЛАУ, так и разные способы хранения вычисляемых коэффициентов.

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

      subroutine prorep (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 Масштабируемость алгоритма и его реализации

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

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

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