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

Встречная прогонка для трёхдиагональной матрицы, точечный вариант
Последовательный алгоритм
Последовательная сложность	[math]8n-2[/math]
Объём входных данных	[math]4n-2[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]2.5n-1[/math]
Ширина ярусно-параллельной формы	[math]4[/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]. Встречная прогонка, как и классическая монотонная, заключается в исключении из уравнений неизвестных, однако, в отличие от монотонной, в ней исключение ведут одновременно с обоих "краёв" СЛАУ (верхнего и нижнего). В принципе, её можно считать простейшим вариантом метода редукции (при m=1 и встречных направлениях монотонных прогонок).

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

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

[math]m[/math] здесь - номер уравнения, на котором "встречаются" две ветви прямого хода - "сверху" и "снизу".

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

"сверху":

[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}), \quad i = 1, 2, \cdots , m-1[/math],

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

и "снизу":

[math]\xi_{N} = a_{N}/c_{N}[/math],

[math]\eta_{N} = f_{N}/c_{N}[/math],

[math]\xi_{i} = a_{i}/(c_{i}-b_{i}\xi_{i+1})[/math], [math]\quad i = N-1, N-2, \cdots , m[/math],

[math]\eta_{i} = (f_{i}+b_{i}\eta_{i+1})/(c_{i}-b_{i}\xi_{i+1})[/math], [math]\quad i = N-1, N-2, \cdots , m. [/math]

после чего вычисляют решение с помощью обратного хода

[math]y_{m} = (\eta_{m}+\xi_{m}\beta_{m})/(1-\xi_{m}\alpha_{m})[/math],

[math]y_{m-1} = (\beta_{m}+\alpha_{m}\eta_{m})/(1-\xi_{m}\alpha_{m})[/math],

[math]y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}, \quad i = m-2, \cdots , 1, 0[/math],

[math]y_{i+1} = \xi_{i+1} y_{i} + \eta_{i+1}, \quad i = m, m+1, \cdots , N-1[/math].

В обычно приводимых формулах встречной прогонки нет формулы для компоненты [math]y_{m-1}[/math], которая вычисляется позже в обратном ходе. Однако это удлиняет критический путь графа, откладывая вычисление [math]y_{m-1}[/math] на момент, когда уже будет вычислена [math]y_{m}[/math], хотя обе компоненты могут быть вычислены одновременно почти независимо друг от друга.

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

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

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

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

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

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

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

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

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

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

[math]\xi_{N} = a_{N}/c_{N}[/math],

[math]\eta_{1} = f_{N}/c_{N}[/math].

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

[math]\alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i})[/math], [math]\quad i = 1, 2, \cdots , m-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 , m-1[/math],

[math]\xi_{i} = a_{i}/(c_{i}-b_{i}\xi_{i+1})[/math], [math]\quad i = N-1, N-2, \cdots , m[/math],

[math]\eta_{i} = (f_{i}+b_{i}\eta_{i+1})/(c_{i}-b_{i}\xi_{i+1})[/math], [math]\quad i = N-1, N-2, \cdots , m[/math].

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

[math]y_{m-1} = (\beta_{m}+\alpha_{m}\eta_{m})/(1-\xi_{m}\alpha_{m})[/math],

[math]y_{m} = (\eta_{m}+\xi_{m}\beta_{m})/(1-\xi_{m}\alpha_{m})[/math].

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

[math]y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}[/math], [math]\quad i = m-1, m-2, \cdots , 1, 0[/math],

[math]y_{i+1} = \xi_{i+1} y_{i} + \eta_{i+1}[/math], [math]\quad i = m, m+1, \cdots , N-1[/math].

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

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

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

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

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

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

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

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

Обе ветви прямого хода можно выполнять одновременно, если [math]N=2m-1[/math], т.е. [math]n=2m[/math]. В этом случае встречная прогонка требует последовательного выполнения следующих ярусов:

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

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

При нечётном [math]n[/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 Свойства алгоритма

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

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

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

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

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

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

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

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

      subroutine vprogm (a,x,N) ! N=2m-1
      real a(3,0:N), x(0:N)

      m=(N+1)/2

      a(2,0)=1./a(2,0)
      a(2,N)=1./a(2,N)

      a(3,0)=-a(3,0)*a(2,0) ! alpha 1
      a(1,N)=-a(1,N)*a(2,N) ! xi N

      x(0)=x(0)*a(2,0) ! beta 1
      x(N)=x(N)*a(2,N) ! eta N

      do 10 i=1,m-1

        a(2,i)=1./(a(2,i)+a(1,i)*a(2,i-1))
        a(2,N-i)=1./(a(2,N-i)+a(3,N-i)*a(2,N-i+1))

        a(3,i) = -a(3,i)*a(2,i) ! alpha i+1
        a(1,N-i) = -a(1,N-i)*a(2,N-i) ! xi N-i

        x(i)=(x(i)-a(1,i)*x(i-1))*a(2,i) ! beta i+1
        x(N-i)=(x(N-i)-a(3,N-i)*x(N-i+1))* a(2,N-i) ! eta N-i

  10  continue

      a(1,0)=1./(1.-a(1,m)*a(3,m-1))
      bb=x(m-1)
      ee=x(m)
      x(m-1)=(bb+a(3,m-1)*ee))*a(1,0) ! y m-1
      x(m) = (ee+a(1,m)*bb)*a(1,0) ! y m

      do 20 i=m+1,N
        x(i)=a(1,i)*x(i-1)+x(i) ! y i
        x(N-i)=a(3,N-i)*x(N-i+1)+x(i) ! y N-i
  20  continue
      return
      end

2.2 Локальность данных и вычислений

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

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

2.2.1 Локальность реализации алгоритма

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

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

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

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

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

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

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

Рисунок 5. Отдельные профили для разных массивов

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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