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

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 Математическое описание алгоритма

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

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

Метод Стоуна в части решения двухдиагональных СЛАУ [math]Ly = b[/math] и [math]D^{-1}Ux = z[/math], полученных при решении исходной [math]Ax = b[/math] после вычисления разложения [math]A = LU[/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} = B_{i} \begin{bmatrix} y_{i-1} \\ 1 \\ \end{bmatrix}, i = 2,..., n [/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] \begin{bmatrix} y_i \\ 1 \\ \end{bmatrix} = B_{i} B_{i-1} ... B_{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_{n} = \begin{bmatrix} s_{n} & t_{n} \\ 0 & 1 \\ \end{bmatrix} = B_{n} B_{n-1} ... B_{2} , 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]z_n = \frac{s_{n}b_{1}+t_{n}}{u_{nn}}[/math],

и окончательные

[math]x_n = z_n[/math],

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

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

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

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

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

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

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

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