Полный метод циклической редукции

Циклическая редукция для трёхдиагональной матрицы, точечный вариант
Последовательный алгоритм
Последовательная сложность	[math]17n + o(n)[/math]
Объём входных данных	[math]4n-2[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(log_2 n)[/math]
Ширина ярусно-параллельной формы	[math]O(n)[/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] A = \begin{bmatrix} b_{1} & a_{1} & 0 & \cdots & \cdots & 0 \\ c_{2} & b_{2} & a_{2} & \cdots & \cdots & 0 \\ 0 & c_{3} & b_{3} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & c_{n-1} & b_{n-1} & a_{n-1} \\ 0 & \cdots & \cdots & 0 & c_{n} & b_{n} \\ \end{bmatrix}\begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \\ \end{bmatrix} = \begin{bmatrix} f_{1} \\ f_{2} \\ \vdots \\ f_{n} \\ \end{bmatrix} [/math]

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

[math]b_{1} x_{1} + a_{1} x_{2} = f_{1}[/math],

[math]c_{i} x_{i-1} + b_{i} x_{i} + a_{i} x_{i+1} = f_{i}, 2 \le i \le n-1[/math],

[math]c_{n} x_{n-1} + b_{n} x_{n} = f_{n}[/math].

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

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

Лучше всего схема циклической редукции разработана для случая [math]n = 2^{m}-1[/math]. Эта схема состоит из прямого и обратного ходов. Прямой ход состоит из последовательного уменьшения в СЛАУ количества уравнений почти в 2 раза (за счёт подстановки из уравнений с нечётными номерами заменяются уравнения с чётными), пока не останется одно уравнение, обратный - в получении всё большего количества компонент решения исходной СЛАУ. Оба хода - как прямой, так и обратный - разбиты на шаги. Здесь мы приведём тот вариант алгоритма, в котором операции экономятся за счёт предварительной нормировки уравнений, используемых для исключения неизвестных.

1.2.1 Прямой ход редукции

В начале считается, что все [math]c^{(0)}_{i} = c_{i}, b^{(0)}_{i} = b_{i}, a^{(0)}_{i} = a_{i}, f^{(0)}_{i} = f_{i}, x^{(0)}_{i} = x_{i}[/math]

На k-м шаге теперь выполняем процедуру редукции системы уравнений размерности n.

Для каждого из уравнений

[math]c^{(k)}_{i} x^{(k)}_{i-1} + b^{(k)}_{i} x^{(k)}_{i} + a^{(k)}_{i} x^{(k)}_{i+1} = f^{(k)}_{i}[/math]

с нечётными [math]i[/math] с помощью деления уравнения на [math]b^{(k)}_{i}[/math] выполняется его замена на уравнение

[math] (c^{(k)}_{i}/b^{(k)}_{i}) x^{(k)}_{i-1} + x^{(k)}_{i} + (a^{(k)}_{i}/b^{(k)}_{i}) x^{(k)}_{i+1} = f^{(k)}_{i}/b^{(k)}_{i}[/math]

Рисунок 1. Микрограф "узла" подготовительного шага прямого хода алгоритма циклической редукции

Уравнение

[math]b^{(k)}_{1} x^{(k)}_{1} + a^{(k)}_{1} x^{(k)}_{2} = f^{(k)}_{1}[/math]

аналогично меняется на уравнение

[math]x^{(k)}_{1} + (a^{(k)}_{1}/b^{(k)}_{1}) x^{(k)}_{2} = f^{(k)}_{1}/b^{(k)}_{1}[/math]:

а уравнение

[math]c^{(k)}_{n} x^{(k)}_{n-1} + b^{(k)}_{n} x^{(k)}_{n} = f^{(k)}_{n}[/math]

меняется на уравнение

[math](c^{(k)}_{n}/b^{(k)}_{n}) x^{(k)}_{n-1} + x^{(k)}_{n} = f^{(k)}_{n}/b^{(k)}_{n}[/math]:

Рисунок 2a. Микрограф "узла" с нечётным номером прямого хода алгоритма циклической редукции

Для каждого же из уравнений

[math]c^{(k)}_{i} x^{(k)}_{i-1} + b^{(k)}_{i} x^{(k)}_{i} + a^{(k)}_{i} x^{(k)}_{i+1} = f^{(k)}_{i}[/math]

с чётными [math]i[/math] (кроме [math]2[/math] и [math]n-2[/math]) выполняется, с учётом [math]x^{(k+1)}_{i/2} = x^{(k)}_{i}[/math] его замена на уравнение

[math]c^{(k+1)}_{i/2} x^{(k+1)}_{(i-2)/2} + b^{(k+1)}_{i/2} x^{(k+1)}_{i/2} + a^{(k+1)}_{i/2} x^{(k+1)}_{(i+2)/2} = f^{(k+1)}_{i/2}[/math]

Рисунок 2b. Микрограф "узла" с чётным номером прямого хода алгоритма циклической редукции

при этом

[math]c^{(k+1)}_{i/2} = - c^{(k)}_{i}(c^{(k)}_{i-1}/b^{(k)}_{i-1})[/math],

[math]a^{(k+1)}_{i/2} = - a^{(k)}_{i}(a^{(k)}_{i+1}/b^{(k)}_{i+1})[/math],

[math]b^{(k+1)}_{i/2} = b^{(k)}_{i} - c^{(k)}_{i}(a^{(k)}_{i-1}/b^{(k)}_{i-1}) - a^{(k)}_{i}(c^{(k)}_{i+1}/b^{(k)}_{i+1})[/math],

[math]f^{(k+1)}_{i/2} = f^{(k)}_{i} - c^{(k)}_{i}f^{(k)}_{i-1}/b^{(k)}_{i-1} - a^{(k)}_{i}f^{(k)}_{i-1}/b^{(k)}_{i-1}[/math].

Для 2го уравнения выполняется его замена на уравнение

[math]b^{(k+1)}_{1} x^{(k+1)}_{1} + a^{(k+1)}_{1} x^{(k+1)}_{(2} = f^{(k+1)}_{1}[/math]

при этом

[math]a^{(k+1)}_{1} = - a^{(k)}_{2}(a^{(k)}_{3}/b^{(k)}_{3})[/math],

[math]b^{(k+1)}_{1} = b^{(k)}_{2} - c^{(k)}_{2}(a^{(k)}_{1}/b^{(k)}_{1}) - a^{(k)}_{2}(c^{(k)}_{3}/b^{(k)}_{3})[/math],

[math]f^{(k+1)}_{1} = f^{(k)}_{2} - c^{(k)}_{2}f^{(k)}_{1}/b^{(k)}_{1} - a^{(k)}_{2}f^{(k)}_{1}/b^{(k)}_{1}[/math]

[math]n-1[/math]-е уравнение заменяется на

[math]c^{(k+1)}_{(n-1)/2} x^{(k+1)}_{(n-3)/2} + b^{(k+1)}_{(n-1)/2} x^{(k+1)}_{(n-1)/2} = f^{(k+1)}_{(n-1)/2}[/math]

при этом

[math]c^{(k+1)}_{(n-1)/2} = - c^{(k)}_{n-1}(c^{(k)}_{n-2}/b^{(k)}_{n-2})[/math],

[math]b^{(k+1)}_{(n-1)/2} = b^{(k)}_{n-1} - c^{(k)}_{n-1}(a^{(k)}_{n-2}/b^{(k)}_{n-2}) - a^{(k)}_{n-1}(c^{(k)}_{n}/b^{(k)}_{n})[/math],

[math]f^{(k+1)}_{(n-1)/2} = f^{(k)}_{n-1} - c^{(k)}_{n-1}f^{(k)}_{n-2}/b^{(k)}_{n-2} - a^{(k)}_{n-1}f^{(k)}_{n-2}/b^{(k)}_{n-2}[/math].

По окончании всех этих манипуляций размерность k+1-й СЛАУ оказывается равной [math](n-1)/2[/math].

Шаги повторяются до тех пор, пока после [math]m-1[/math] шагов редукции размерность СЛАУ не становится равной 1 и остаётся одно уравнение

[math]b^{(m-1)}_{1} x^{(m-1)}_{1} = f^{(m-1)}_{1}[/math]

1.2.2 Обратный ход редукции

Рисунок 3. Микрограф "узла" обратного хода алгоритма циклической редукции

Из последнего уравнения, полученного прямым ходом, вычисляется

[math]x^{(m-1)}_{1} = f^{(m-1)}_{1}/b^{(m-1)}_{1}[/math]

Теперь, последовательно уменьшая верхние индексы неизвестных, используется нечётные уравнения каждого шага для вычисления неизвестных с соотвествующими нечётными номерами. Чётные неизвестные получаются из тождеств [math]x^{(k)}_{i} = x^{(k+1)}_{i/2}[/math], а для нечётных [math]i[/math]

[math]x^{(k)}_{i} = f^{(k)}_{i}/b^{(k)}_{i} - (c^{(k)}_{i}/b^{(k)}_{i}) x^{(k)}_{i-1} - (a^{(k)}_{i}/b^{(k)}_{i}) x^{(k)}_{i+1} [/math]

c "левого края" системы будет

[math]x^{(k)}_{1} = f^{(k)}_{1}/b^{(k)}_{1} - (a^{(k)}_{1}/b^{(k)}_{1}) x^{(k)}_{2}[/math]

а с "правого"

[math]x^{(k)}_{n} = f^{(k)}_{n}/b^{(k)}_{n} - (c^{(k)}_{n}/b^{(k)}_{n}) x^{(k)}_{n-1}[/math]

После вычисления всех [math]x^{(0)}_{i}[/math] значения искомых неизвестных [math]x_{i} = x^{(k)}_{i}[/math].

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

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

[math] c^{(k)}_{i}/b^{(k)}_{i} , a^{(k)}_{i}/b^{(k)}_{i} , f^{(k)}_{i}/b^{(k)}_{i}[/math]

почти все используются для преобразований двух чётных уравнений, то при выделении "микровычислений", из которых следует составить шаги редукции и которые составляют его ядро, лучше отнести вычисления этих отношений к предыдущему шагу редукции. Таким образом, на "подготовительном шаге" микроядро будет для каждого нечётного [math]i[/math] состоять только из трёх делений (кроме [math]2[/math] и [math]n[/math] - там будет по 2 деления), а затем на каждом последующем шаге редукции для каждого нечётного [math]i[/math] - из двух умножений и двух вычислений выражений типа [math]a-bc-de[/math], с последующими тремя делениями (на "краях" часть этих операций отсутствует или урезана, но общую картину это не очень меняет). Для чётных [math]i[/math] деления в микроядре будут отсутствовать.

Что касается шагов обратного хода, то там для каждого [math]i[/math] рано или поздно выполняется одна операция типа [math]a-bc-de[/math] (на "краях" - типа [math]a-bc[/math]).

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

Рисунок 4. Граф алгоритма прямого хода циклической редукции при n=15. Светло-зелёным обозначены микроядра "подготовительного шага" с тремя (или меньше) делениями, синим - микроядра нечётных узов, сине-зелёным - ядра чётных узлов (без делений).

Рисунок 5. Граф алгоритма обратного хода циклической редукции при n=15. В вершинах - операции вида a-bc-de, в вершинах с 1 входящей дугой - вида a-bc. Показаны только дуги, передающие значения найденных неизвестных.

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

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

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

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

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

Поскольку алгоритм циклической редукции обычно не предназначен для последовательного исполнения, то его "последовательную сложность" скорее следует считать только оценкой общего количества операций. Если взять только главный член, линейный по размеру, и опустить меньшие по порядку (логарифмические и константы), то получается, что в циклической редукции 3n делений, 8n умножений и 6n вычитаний/сложений. Таким образом, при классификации по последовательной сложности, алгоритм циклической редукции относится к алгоритмам с линейной сложностью.

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

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

У прямого хода циклической редукции макрограф представлен на Рис. 4, при этом графы разные макровершин представлены на Рис. 1, 2a, 2b. Как видно, после первого шага схему можно условно назвать "схемой страивания", поскольку связанные с каждым уравнением, остающимся в редуцированной СЛАУ, коэффициенты вычисляются по формулам из коэффициентов трёх соседних уравнений в нередуцированной СЛАУ.

Обратный ход имеет макрограф, представленный на Рис. 5, при этом графы макровершин представлены на Рис. 3 (крайние макровершины, имеющие только одну входящую дугу на макрографе, замещают недостающие данные нулями, соответственно уменьшая вычисления). Этот макрограф, в принципе, имеет макроструктуру, обратную к схеме сдваивания, поэтому её можно условно назвать "схемой раздваивания", с одной поправкой: все макровершины, кроме крайних, имеют по две входящих дуги, которые немного изменяют обычное "раздваивание".

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

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

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

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

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

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

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

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

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

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