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

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

Встречная повторная прогонка, как и повторная монотонная, заключается в исключении из уравнений неизвестных, однако, в отличие от монотонной, в ней исключение ведут одновременно с обоих "краёв" СЛАУ (верхнего и нижнего). В принципе, её можно считать простейшим вариантом метода редукции (при m=1 и встречных направлениях монотонных прогонок).

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

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

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

"сверху":

[math]\beta_{1} = f_{0}/c_{0}, [/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]\eta_{N} = f_{N}/c_{N}[/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]1/c_{0}[/math], [math]1/c_{N}[/math], [math]1/(c_{i}-a_{i}\alpha_{i})[/math], [math]1/(c_{i}-b_{i}\xi_{i+1})[/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], хотя обе компоненты могут быть вычислены одновременно почти независимо друг от друга, так и в случае с нечётным числом переменных, когда для вычисления на "месте встречи" нужно подождать дополнительно одно вычисление коэффициентов либо "сверху" от него, либо "снизу".

Поэтому в более поздних источниках приводятся формулы, которые более оптимальны для нечётного количества неизвестных. В них старт обратного хода заменяется на формулу (в наших обозначениях)

[math]y_{m} = (f_{m}+b_{m}\eta_{m+1}+a_{m}\beta_{m})/(c_{m}-a_{m}\alpha_{m}-b_{m}\xi_{m+1})[/math]

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

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

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

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

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

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

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

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

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

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

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

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

[math]\eta_{i} = (f_{i}+b_{i}\eta_{i+1})(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/(1-\xi_{m}\alpha_{m}))[/math],

[math]y_{m} = (\eta_{m}+\xi_{m}\beta_{m})(1/(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].

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

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

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

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

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

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

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

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

• по [math]3m-2[/math] ярусов умножений и [math]2m-1[/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).

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

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

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

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

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

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

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

      subroutine vprogmr (a,x,N) ! N=2m-1
      real a(3,N), x(N)

      m=(N+1)/2

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

      do 10 i=2,m-1

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

  10  continue

      x(m) =(x(m)-a(1,m)*x(m-1)-a(3,m)*x(m+1))*a(2,m) ! y m

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

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

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

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

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

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

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

Из рисунка можно предположить, что переборы, скорее всего, устроены достаточно регулярно и не меняют своего поведения в течение программы. Однако нужно детально разобрать, как устроены эти переборы; для этого рассмотрим некоторые наиболее интересные фрагменты профиля более детально.

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

Рисунок 3. Профиль обращений, фрагмент 1

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

Рисунок 4. Профиль обращений, фрагмент 2

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

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

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

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

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

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

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

С появлением суперскалярных процессоров оказалось, что для получения выигрыша (около 20%) перед монотонной повторной прогонкой встречную даже необязательно "раскидывать" на 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.4 Масштабируемость алгоритма и его реализации

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

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

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

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

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

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

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