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

Прогонка для трёхдиагональной матрицы, точечный вариант
Последовательный алгоритм
Последовательная сложность	$8n-2$
Объём входных данных	$4n-2$
Объём выходных данных	$n$
Параллельный алгоритм
Высота ярусно-параллельной формы	$2.5n+3$
Ширина ярусно-параллельной формы	$4$

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

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

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

Встречная прогонка - один из вариантов метода исключения неизвестных в приложении к решению СЛАУ вида $Ax = b$, где

$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}$

Часто, однако, при изложении сути метода прогонки элементы правой части и матрицы системы обозначают и нумеруют по-другому, например СЛАУ может иметь вид ($N=n-1$)

$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}$

или, если записывать отдельно по уравнениям, то

$c_{0} y_{0} - b_{0} y_{1} = f_{0},\\ -a_{i} y_{i-1} + c_{i} y_{i} - b_{i} y_{i+1} = f_{i}, 1 \le i \le N-1, \\ -a_{N} y_{N-1} + c_{N} y_{N} = f_{N}$

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

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

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

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

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

"сверху":

$\alpha_{1} = b_{0}/c_{0},\\ \beta_{1} = f_{0}/c_{0}, \\ \alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i}), \quad i = 1, 2, \cdots , m-1, \\ \beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i}), \quad i = 1, 2, \cdots , m-1.$

и "снизу":

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

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

$y_{m} = (\eta_{m}+\xi_{m}\beta_{m})/(1-\xi_{m}\alpha_{m}), \\ y_{m-1} = (\beta_{m}+\alpha_{m}\eta_{m})/(1-\xi_{m}\alpha_{m}), \\ y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}, \quad i = m-2, \cdots , 1, 0, \\ y_{i+1} = \xi_{i+1} y_{i} + \eta_{i+1}, \quad i = m, m+1, \cdots , N-1.$

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

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

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

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

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

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

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

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

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

$\alpha_{1} = b_{0}/c_{0},\\ \beta_{1} = f_{0}/c_{0}, \\ \xi_{N} = a_{N}/c_{N},\\ \eta_{1} = f_{N}/c_{N}.$

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

$\alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i}), \quad i = 1, 2, \cdots , m-1, \\ \beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i}), \quad i = 1, 2, \cdots , m-1, \\ \xi_{i} = a_{i}/(c_{i}-b_{i}\xi_{i+1}), \quad i = N-1, N-2, \cdots , m, \\ \eta_{i} = (f_{i}+b_{i}\eta_{i+1})/(c_{i}-b_{i}\xi_{i+1}), \quad i = N-1, N-2, \cdots , m.$

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

$y_{m-1} = (\beta_{m}+\alpha_{m}\eta_{m})/(1-\xi_{m}\alpha_{m}), \\ y_{m} = (\eta_{m}+\xi_{m}\beta_{m})/(1-\xi_{m}\alpha_{m}).$

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

$y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}, \quad i = m-1, m-2, \cdots , 1, 0, \\ y_{i+1} = \xi_{i+1} y_{i} + \eta_{i+1}, \quad i = m, m+1, \cdots , N-1.$

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

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

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

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

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

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

Обе ветви прямого хода выполняются одновременно для $N=2m-1$, т.е. для $n=2m$. Тогда для выполнения встречной прогонки требуется последовательно выполнить следующие ярусы:

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

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

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

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

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

Объём выходных данных: $n$.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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