Алгоритм сдваивания Стоуна для 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} \\ 1 & 0 \\ \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.2 Локальность данных и вычислений

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

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

Из-за неустойчивости алгоритма любые способы реализации теряют смысл.

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

Из-за неустойчивости алгоритма не имеют смысла любые разговоры о масштабируемости любых его реализаций.

2.5 Динамические характеристики и эффективность реализации алгоритма

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

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

На области относительно приемлемых ошибок (малые размеры задач) у метода всё равно нет шансов на нормальную реализацию - быстрее будет СЛАУ прогонкой на обычной персоналке решить.

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

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