Метод сдваивания Стоуна для решения двудиагональных СЛАУ

Метод сдваивания Стоуна для решения двудиагональных СЛАУ
Последовательный алгоритм
Последовательная сложность	[math]3(n-1)(\lceil \log_2 (n-1) \rceil+2)[/math]
Объём входных данных	[math]4n-2[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]4 \lceil \log_2 (n-1) \rceil + 5[/math]
Ширина ярусно-параллельной формы	[math]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]

Метод сдваивания Стоуна впервые предложен в начале 70-х гг. 20го века в качестве альтернативы другим параллельным алгоритмам решения трёхдиагональных СЛАУ, например, методу циклической редукции.

Здесь рассматривается его вторая часть - решение двух двудиагональных СЛАУ. Оно использует представление матрицы

[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} [/math]

в виде произведения матриц

[math] L = \begin{bmatrix} 1 & 0 & 0 & \cdots & \cdots & 0 \\ l_{21} & 1 & 0 & \cdots & \cdots & 0 \\ 0 & l_{32} & 1 & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & l_{n-1 n-2} & 1 & 0 \\ 0 & \cdots & \cdots & 0 & l_{n n-1} & 1 \\ \end{bmatrix} [/math]

и

[math] U = \begin{bmatrix} u_{11} & u_{12} & 0 & \cdots & \cdots & 0 \\ 0 & u_{22} & u_{23}& \cdots & \cdots & 0 \\ 0 & 0 & u_{33} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & 0 & u_{n-1 n-1} & u_{n-1 n} \\ 0 & \cdots & \cdots & 0 & 0 & u_{n n} \\ \end{bmatrix} [/math]

Важным моментом является то, что алгоритм Стоуна использует то же самое разложение, что вычисляется не только в первой части метода (алгоритме сдваивания Стоуна для LU-разложения трёхдиагональной матрицы), но и в устойчивой компактной схеме метода Гаусса.

При уже полученном разложении матрицы решение СЛАУ [math]Ax = b[/math] можно поменять на последовательное решение двух СЛАУ [math]Ly = b[/math] и затем [math]Ux = y[/math]. При этом вторую СЛАУ тоже можно решить как последовательность СЛАУ [math]Dz = y[/math] и [math]D^{-1}Ux = z[/math], где [math]D[/math] - диагональная матрица, составленная из диагональных элементов матрицы [math]U[/math].

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

Метод Стоуна в части решения двухдиагональных СЛАУ [math]Ly = b[/math] и [math]D^{-1}Ux = z[/math], полученных при решении исходной [math]Ax = b[/math] после вычисления разложения [math]A = LU[/math], заключается в том, что получающиеся при их непосредственном решении рекурсивные зависимости

Рисунок 1. Граф вычисления матриц [math]K_{i}[/math] при [math]n=9[/math]. Каждая вершина соответствует составной операции, состоящей из одного умножения и одной операции [math]a+bc[/math]. В чёрных вершинах вычисляются используемые результаты ([math]s_{i}[/math] и [math]t_{i}[/math]), в светлых - промежуточные

Рисунок 2. Граф вычисления матриц [math]R_{i}[/math] при [math]n=9[/math]. Каждая вершина соответствует составной операции, состоящей из одного умножения и одной операции [math]a+bc[/math]. В чёрных вершинах вычисляются используемые результаты ([math]v_{i}[/math] и [math]w_{i}[/math]), в светлых - промежуточные

[math]y_1 = b_1[/math], и потом [math]y_{i} = b_{i} - l_{i i-1} y_{i-1}, i = 2,..., n[/math]

и

[math]x_n = z_n[/math], и потом [math]x_{i} = z_{i} - \frac{u_{i i+1}}{u_{ii}} x_{i+1}, i = n-1,...,1[/math]

заменяются соответственно на

[math] \begin{bmatrix} y_i \\ 1 \\ \end{bmatrix} = \begin{bmatrix} l_{i i-1} & b_{i} \\ 0 & 1 \\ \end{bmatrix} \begin{bmatrix} y_{i-1} \\ 1 \\ \end{bmatrix} = P_{i} \begin{bmatrix} y_{i-1} \\ 1 \\ \end{bmatrix}, i = 2,..., n [/math], где [math] P_{i} = \begin{bmatrix} l_{i i-1} & b_{i} \\ 0 & 1 \\ \end{bmatrix} [/math]

и

[math] \begin{bmatrix} x_i \\ 1 \\ \end{bmatrix} = \begin{bmatrix} \frac{u_{i i+1}}{u_{ii}} & z_{i} \\ 0 & 1 \\ \end{bmatrix} \begin{bmatrix} x_{i+1} \\ 1 \\ \end{bmatrix} = C_{i} \begin{bmatrix} x_{i+1} \\ 1 \\ \end{bmatrix}, i = n-1,...,1[/math], где [math] C_{i} = \begin{bmatrix} \frac{u_{i i+1}}{u_{ii}} & z_{i} \\ 0 & 1 \\ \end{bmatrix} [/math]

и после выполнения подстановок оказывается, что

[math] \begin{bmatrix} y_i \\ 1 \\ \end{bmatrix} = P_{i} P_{i-1} ... P_{2} \begin{bmatrix} b_{1} \\ 1 \\ \end{bmatrix}, i = 2,..., n [/math],

[math] \begin{bmatrix} x_i \\ 1 \\ \end{bmatrix} = C_{i} C_{i+1} ... C_{n-1} \begin{bmatrix} z_{n} \\ 1 \\ \end{bmatrix}, i = 1,...,n-1[/math]

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

[math] K_{i} = \begin{bmatrix} s_{i} & t_{i} \\ 0 & 1 \\ \end{bmatrix} = P_{i} P_{n-1} ... P_{2} , i = 2,..., n [/math]

После этого вычисляются промежуточные результаты

[math]z_i = \frac{s_{i}b_{1}+t_{i}}{u_{ii}}, i = 2,..., n[/math],

и снова по схеме сдваивания матрицы

[math] R_{i} = \begin{bmatrix} v_{i} & w_{i} \\ 0 & 1 \\ \end{bmatrix} = C_{i} C_{i+1} ... C_{n-1} , i = 1,...,n-1 [/math].

После чего окончательно вычисляются

[math]x_n = z_n[/math],

[math]x_i = v_{i}z_{n}+w_{i}, i = 1,...,n-1 [/math].

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

Основную часть вычислительного ядра алгоритма составляют операции трёх типов: деление, умножение и операция типа [math]a+bc[/math].

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

На макроуровне можно выделить такие 4 макрооперации: вычисление сдваиванием матриц [math]K_{i}[/math], вычисление элементов матриц [math]C_i[/math], вычисление сдваиванием матриц [math]R_{i}[/math], вычисление результатов.

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

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

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

Для полного выполнения алгоритма Стоуна и решения двух двухдиагональных СЛАУ нужно выполнить (если не используются предвычисления обратных элементов):

• [math]2n-2[/math] делений,
• [math]2(n-1)(\lceil \log_2 (n-1) \rceil+1)[/math] умножений,
• [math](n-1)(\lceil \log_2 (n-1) \rceil+2)[/math] сложений.

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

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

Как уже отмечено, макроструктура алгоритма состоит из 4 частей.

В первой части производится вычисление сдваиванием матриц [math]K_{i}[/math]. Граф этой части показан на Рисунке 1.

Во второй части происходит вычисление элементов матриц [math]C_i[/math]. При этом внедиагональные элементы ([math]z_i[/math]) вычисляются с использованием результатов первой части, а вот диагональные (частные [math]\frac{u_{i i+1}}{u_{ii}}[/math]) можно вычислить заранее, одновременно с 1й частью.

В третьей части - вычисление сдваиванием матриц [math]R_{i}[/math], оно показано на рисунке 2 и использует результаты 2й части.

В четвёртой, последней, части - вычисление результатов, с использованием результатов 2й части.

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

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

На критическом пути алгоритма Стоуна решения двух двухдиагональных СЛАУ нужно выполнить (если не используются предвычисления обратных элементов):

• [math]1[/math] ярус делений,
• [math]2(\lceil \log_2 (n-1) \rceil+2)[/math] ярусов умножений,
• [math]2(\lceil \log_2 (n-1) \rceil+2)[/math] ярусов сложений.

Поэтому алгоритм должен быть отнесён к алгоритмам логарифмической сложности по количеству последовательных операций. Ширина яруса равна [math]n[/math], поэтому алгоритм должен быть отнесён к алгоритмам линейной сложности по ширине ярусов.

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

Входные данные: двудиагональные матрицы [math]L[/math] (с единичной диагональю) (элементы [math]l_{ij}[/math]) и двудиагональные матрицы [math]U[/math] (элементы [math]u_{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.2 Локальность данных и вычислений

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

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

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

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

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