Алгоритм сдваивания Стоуна для LU-разложения трёхдиагональной матрицы

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

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

Алгоритм сдваивания Стоуна для LU-разложения трёхдиагональной матрицы - часть метода сдваивания Стоуна для решения СЛАУ вида [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]LU[/math]-разложение. Оно состоит в представлении матрицы

Рисунок 1. Граф вычисления произведений матриц [math]T_{i}[/math] при [math]n=9[/math]. Каждая вершина на первом ярусе соответствует составной операции, состоящей из одного умножения и одной операции [math]a+bc[/math], на остальных ярусах - из двух операций [math]ab+cd[/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} [/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]

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

Теоретически метод Стоуна основан на том, что, если ввести величины [math]q_i[/math] так, что [math]q_0 = 1[/math], [math]q_1 = a_{11}[/math] и [math]q_i = a_{ii} q_{i-1} - a_{ii-1}a_{i-1i} q_{i-2}[/math], то после вычисления всех [math]q_i[/math] легко вычислить все [math]u_{ii} = q_i/q_{i-1}[/math], [math]l_{ii-1} = a_{ii-1}/u_{i-1i-1}[/math]. При этом не нужно вычислять [math]u_{ii+1} = a_{ii+1}[/math].

Вычисление всех [math]q_i[/math] производится с использованием ассоциативности матричного умножения. Из рекуррентной связи [math]q_i = a_{ii} q_{i-1} - a_{ii-1}a_{i-1i} q_{i-2}[/math] следует матричное равенство

[math] \begin{bmatrix} q_i \\ q_{i-1} \\ \end{bmatrix} = \begin{bmatrix} a_{ii} & -a_{ii-1}a_{i-1i} \\ 0 & 1 \\ \end{bmatrix} \begin{bmatrix} q_{i-1} \\ q_{i-2} \\ \end{bmatrix} = T_i \begin{bmatrix} q_{i-1} \\ q_{i-2} \\ \end{bmatrix} = T_i T_{i-1}...T_2 \begin{bmatrix} q_{1} \\ q_{0} \\ \end{bmatrix} = T_i T_{i-1}...T_2 \begin{bmatrix} a_{11} \\ 1 \\ \end{bmatrix} [/math]

Произведения матриц [math]T_i T_{i-1}...T_2[/math] вычисляются параллельной схемой сдваивания. При этом благодаря тому, что вторая строка всех матриц [math]T_i[/math] равна строке единичной матрицы, вторая строка произведения [math]T_i T_{i-1}...T_k[/math] совпадает с первой строкой произведения [math]T_{i-1}...T_k[/math], что позволяет сэкономить вычисления вдвое против перемножения неособенных матриц этого же порядка.

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

Вначале полагается для всех [math]i[/math] от 2 до [math]n[/math]

[math]p_i = a_{ii}, r_i = -a_{ii-1}a_{i-1i}[/math],

а также

[math]p_1 = 0, r_1 = 1[/math].

Затем вычисления на каждом шаге (номер шага обозначим [math]k[/math]) ведутся так:

начиная с [math]i[/math] с 2, пропускают [math]2^{k-1}[/math] чисел, а потом для [math]2^{k-1}[/math] чисел берут ближайшее для них снизу [math]j[/math] вида [math]j = 1 + 2^{k} m [/math].

Для первого шага перевычисляются новые значения [math]p_i, r_i[/math] по формулам: [math]p^{new}_i = p_i p_{i-1}, r^{new}_i = p_i r_{i-1} + r_i[/math], а для остальных шагов - по формулам [math]p^{new}_i = p_i p_{j}+r_i p_{j-1}, r^{new}_i = p_i r_{j} + r_i r_{j-1}[/math].

Затем опять пропускают [math]2^{k-1}[/math] чисел, повторяют для [math]2^{k-1}[/math] чисел такие же вычисления, и так до исчерпания диапазона всего значений [math]i[/math] от 2 до [math]n[/math].

Шаги повторяют до тех пор, пока не оказывается, что нужно пропустить все значения [math]i[/math] от 2 до [math]n[/math]. После этого выполняют вычисления

[math]q_i = p_i a_{11} + r_i[/math] для всех [math]i[/math] от 1 до [math]n[/math] и затем [math]u_{ii} = q_i/q_{i-1}[/math], [math]l_{ii-1} = a_{ii-1}/u_{i-1i-1}[/math] для всех [math]i[/math] от 2 до [math]n[/math].

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

Вычислительное ядро алгоритма составляют в основной части операции типа [math]ab+cd[/math] и [math]ab+c[/math], с небольшой добавкой отдельных умножений и делений. Изолированные длинные ветви вычислений отсутствуют.

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

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

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

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

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

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

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