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

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

Основные авторы описания: А.В.Фролов, Вад.В.Воеводин (раздел 2.2), А.М.Теплов (раздел 2.4)

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=8 без отображения входных и выходных данных. / - деление, f - операция a+bc или a-bc.

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

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

[math]\alpha_{1} = b_{0}/c_{0}[/math],

[math]\beta_{1} = f_{0}/c_{0}[/math],

[math]\alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i})[/math], [math]\quad i = 1, 2, \cdots , N-1[/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 Вычислительное ядро алгоритма

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

Рисунок 2. Детальный граф алгоритма прогонки с однократным вычислением обратных чисел при n=4 без отображения входных и выходных данных. inv - вычисление обратного числа, mult - операция перемножения чисел. Серым цветом выделены операции, повторяющиеся при замене правой части СЛАУ

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

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

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

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

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

[math]\alpha_{1} = b_{0}/c_{0}[/math],

[math]\beta_{1} = f_{0}/c_{0} [/math]

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

[math]\alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i})[/math],

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

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

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

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

• [math]2n-1[/math] делений,
• [math]3n-3[/math] сложений/вычитаний,
• [math]3n-3[/math] умножений.

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

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

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

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

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

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

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

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

Входные данные: трёхдиагональная матрица [math]A[/math] (элементы [math]a_{ij}[/math]), вектор [math]b[/math] (элементы [math]b_{i}[/math]).

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

Объём выходных данных: [math]n[/math].

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

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

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

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

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

В случае, когда требуется решение нескольких СЛАУ с одной и той же матрицей, левую ветвь вычислений (см. рисунки с графом алгоритма) можно не повторять. Это связано с тем, что [math]LU[/math]-разложение матрицы системы не нужно перевычислять. Тогда будет более предпочтителен вариант с заменой делений.

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

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

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

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

      subroutine progm (a,x,N)
      real a(3,0:N), x(0:N)
      a(2,0)=1./a(2,0)
      a(3,0)=-a(3,0)*a(2,0) ! alpha 1
      x(0)=x(0)*a(2,0) ! beta 1
      do 10 i=1,N-1
        a(2,i)=1./(a(2,i)+a(1,i)*a(2,i-1))
        a(3,i)=-a(3,i)*a(2,i) ! alpha i+1
        x(i)=(x(i)-a(1,i)*x(i-1))*a(2,i) ! beta i+1
  10  continue
      a(2,N)=1./(a(2,N)+a(1,N)*a(2,N-1))
      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.2.1 Локальность реализации алгоритма

2.2.1.1 Структура обращений в память и качественная оценка локальности

Рисунок 3. Прогонка, точечный вариант. Общий профиль обращений в память

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

Рисунок 4. Фрагмент общего профиля (выделенный зеленым цветом на рис. 3)

Однако при более детальном рассмотрении можно заметить, что профиль устроен несколько сложнее (рис. 4). Для некоторых массивов на каждом шаге выполняется по несколько очень близких обращений, что приводит к улучшению как пространственной, так и временной локальности.

Таким образом, общий профиль обладает высокой пространственной локальностью (поскольку здесь доминируют последовательные переборы элементов) и средней временной (поскольку некоторые элементы сразу используются повторно).

2.2.1.2 Количественная оценка локальности

Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен здесь (функция Kernel). Условия запуска описаны здесь.

Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.

На рисунке 5 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Для данной программы значения daps является достаточно неожиданным. Согласно качественному анализу, локальность обращений в память в этом случае достаточно высока. Однако daps оценивает производительность работы с памятью значение ниже, чем для теста RandomAccess! Такое различие между локальностью и производительностью встречается нечасто.

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

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

Рисунок 5. Сравнение значений оценки daps

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

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

На рисунке 6 значения приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). В данном случае все предыдущие выводы подтверждаются – локальность оказывается достаточно высокой, что и было замечено согласно качественной оценке локальности.

Рисунок 6. Сравнение значений оценки cvg

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

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

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

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

Проведём исследование масштабируемости реализации алгоритма согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета.

Набор и границы значений изменяемых параметров запуска реализации алгоритма:

• число процессоров 1;
• размер области [10240 : 1024000] с шагом 10240.

В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:

• минимальная эффективность реализации 0.019%;
• максимальная эффективность реализации 0.0255%.

На следующих рисунках приведены графики производительности и эффективности выбранной реализации алгоритма в зависимости от изменяемых параметров запуска.

Рисунок 7. Реализация алгоритма. Изменение производительности в зависимости от размера вектора.

Рисунок 8. Реализация алгоритма. Изменение эффективности в зависимости от размера вектора.

Малая эффективность, по-видимому, связана с избыточной локальностью, описанной в разделе о локальности данных и вычислений.

Исследованная реализация алгоритма

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

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

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

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

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

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