Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][досмотренная версия]
 
(не показано 47 промежуточных версий 3 участников)
Строка 1: Строка 1:
 +
{{algorithm
 +
| name              = Алгоритм сдваивания Стоуна для LU-разложения трёхдиагональной матрицы
 +
| serial_complexity = <math>3(n-1)\lceil \log_2 (n-1) \rceil + o((n-1)\lceil \log_2 (n-1) \rceil)</math>
 +
| pf_height        = <math>2 \lceil \log_2 (n-1) \rceil + 3</math>
 +
| pf_width          = <math>n</math>
 +
| input_data        = <math>3n-2</math>
 +
| output_data      = <math>2n-1</math>
 +
}}
 +
 +
Основные авторы описания: [[Участник:Frolov|А.В.Фролов]]
 +
 
== Свойства и структура алгоритмов ==
 
== Свойства и структура алгоритмов ==
 
=== Общее описание алгоритма ===
 
=== Общее описание алгоритма ===
  
 
'''Алгоритм сдваивания Стоуна для LU-разложения трёхдиагональной матрицы''' - часть [[Метод сдваивания Стоуна|метода сдваивания Стоуна]] для решения СЛАУ<ref name="VOLA">Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.</ref><ref name="MIV">Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.</ref> вида <math>Ax = b</math>, где  
 
'''Алгоритм сдваивания Стоуна для LU-разложения трёхдиагональной матрицы''' - часть [[Метод сдваивания Стоуна|метода сдваивания Стоуна]] для решения СЛАУ<ref name="VOLA">Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.</ref><ref name="MIV">Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.</ref> вида <math>Ax = b</math>, где  
 +
 +
{{Шаблон:Трёхдиагональная СЛАУ в стандартном виде}}
 +
 +
[[Метод сдваивания Стоуна]] впервые предложен в начале 70-х гг. 20го века<ref name="STONE">Stone H.S. An Efficient Parallel Algorithm for the Solution of a Tridiagonal Linear System of Equations // J. ACM, Vol. 20, No. 1 (Jan. 1973), P. 27-38.</ref> в качестве альтернативы другим параллельным алгоритмам решения трёхдиагональных СЛАУ, например, [[Полный метод циклической редукции|методу циклической редукции]].
 +
 +
Здесь рассматривается его первая часть - <math>LU</math>-разложение. Оно состоит в представлении матрицы
 +
 +
{{Шаблон:Трёхдиагональная матрица стандартной формы}}
 +
 +
в виде произведения матриц
 +
 +
{{Шаблон:L2dUniDiag}}
 +
 +
и
 +
 +
{{Шаблон:U2dCommon}}
 +
 +
Важным моментом является то, что в условиях точных вычислений алгоритм сдваивания Стоуна вычисляет то же самое разложение, что и [[Компактная схема метода Гаусса для трёхдиагональной матрицы, последовательный вариант|компактная схема метода Гаусса]].
 +
 +
Теоретически метод Стоуна основан на том, что, если ввести величины <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>
 
:<math>
−
A = \begin{bmatrix}
+
\begin{bmatrix}
−
  a_{11} & a_{12}  & 0 &    \cdots & \cdots & 0 \\
+
  q_i \\
−
a_{21} & a_{22} & a_{23}& \cdots & \cdots & 0 \\
+
  q_{i-1} \\
−
0 &  a_{32} & a_{33} &    \cdots & \cdots & 0 \\
+
\end{bmatrix} = \begin{bmatrix}
−
  \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\
+
  a_{ii}  & -a_{ii-1}a_{i-1i} \\
−
0 & \cdots & \cdots & a_{n-1 n-2} & a_{n-1 n-1} & a_{n-1 n} \\
+
  1 & 0 \\
−
  0 & \cdots & \cdots & 0 & a_{n n-1}  & a_{n n} \\
+
\end{bmatrix} \begin{bmatrix}
−
\end{bmatrix}, x = \begin{bmatrix}
+
  q_{i-1} \\
−
x_{1} \\
+
  q_{i-2} \\
−
x_{2} \\
+
\end{bmatrix} = T_i \begin{bmatrix}
−
\vdots \\
+
  q_{i-1} \\
−
x_{n} \\
+
q_{i-2} \\
−
\end{bmatrix}, b = \begin{bmatrix}
+
\end{bmatrix} = T_i T_{i-1}...T_2 \begin{bmatrix}
−
b_{1} \\
+
  q_{1} \\
−
b_{2} \\
+
q_{0} \\
−
\vdots \\
+
\end{bmatrix} = T_i T_{i-1}...T_2 \begin{bmatrix}
−
b_{n} \\
+
  a_{11} \\
 +
1 \\
 
\end{bmatrix}
 
\end{bmatrix}
 
</math>
 
</math>
  
−
[[Метод сдваивания Стоуна]] впервые предложен в начале 70-х гг. 20го века<ref>Stone H.S. An Efficient Parallel Algorithm for the Solution of a Tridiagonal Linear System of Equations // J. ACM, Vol. 20, No. 1 (Jan. 1973), P. 27-38.</ref> в качестве альтернативы другим параллельным алгоритмам решения трёхдиагональных СЛАУ, например, [[Метод циклической редукции|методу циклической редукции]].
+
Произведения матриц <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>, что позволяет сэкономить вычисления вдвое против перемножения неособенных матриц этого же порядка.
 +
 
 +
[[file:DoublingPartUniver.png|thumb|right|600px|Рисунок 1. Граф вычисления произведений матриц <math>T_{i}</math> при <math>n=9</math>. Каждая вершина на первом ярусе соответствует составной операции, состоящей из одного умножения и одной операции <math>a+bc</math>, на остальных ярусах - из двух операций <math>ab+cd</math>. В чёрных вершинах вычисляются используемые далее результаты, в светлых - промежуточные]]
 +
 
 +
=== Математическое описание алгоритма ===
  
−
Здесь рассматривается его первая часть - <math>LU</math>-разложение.  
+
Вначале полагается для всех <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>.
  
−
=== Математическое описание ===
 
 
=== Вычислительное ядро алгоритма ===
 
=== Вычислительное ядро алгоритма ===
 +
 +
Вычислительное ядро алгоритма составляют в основной части операции типа <math>ab+cd</math> и <math>ab+c</math>, с небольшой добавкой отдельных умножений и делений. Изолированные длинные ветви вычислений отсутствуют.
 +
 
=== Макроструктура алгоритма ===
 
=== Макроструктура алгоритма ===
−
=== Описание схемы реализации последовательного алгоритма ===
+
 
 +
На макроуровне можно выделить такие 4 макрооперации: вычисление элементов матриц <math>T_i</math>, вычисление сдваиванием произведений матриц <math>T_{i}</math>, вычисление результатов.
 +
 
 +
=== Схема реализации последовательного алгоритма ===
 +
 
 +
Метод Стоуна изначально спроектирован для параллельного исполнения, поскольку является по отношению к, например, классической прогонке, алгоритмом с избыточными вычислениями. Смысла в его последовательной реализации нет ещё и из-за того, что он неустойчив.
 +
 
 
=== Последовательная сложность алгоритма ===
 
=== Последовательная сложность алгоритма ===
 +
 +
Для полного выполнения алгоритма Стоуна и получения разложения трёхдиагональной матрицы на две двудиагональные нужно выполнить:
 +
 +
* <math>2n-2</math> делений,
 +
* <math>2(n-1)\lceil \log_2 (n-1) \rceil + o((n-1)\lceil \log_2 (n-1) \rceil)</math> умножений,
 +
* <math>(n-1)\lceil \log_2 (n-1) \rceil +o((n-1)\lceil \log_2 (n-1) \rceil)</math> сложений.
 +
 +
Поэтому алгоритм должен быть отнесён к алгоритмам ''линейно-логарифмической сложности'' по количеству последовательных операций.
 +
 
=== Информационный граф ===
 
=== Информационный граф ===
−
=== Описание ресурса параллелизма алгоритма ===
 
−
=== Описание входных и выходных данных ===
 
−
=== Свойства алгоритма===
 
  
−
== Программная реализация алгоритмов ==
+
Как уже отмечено, макроструктура алгоритма состоит из 4 частей.
 +
 
 +
В первой части производится вычисление недиагональных элементов матриц <math>T_{i}</math>. 
 +
 
 +
Во второй части - вычисление сдваиванием произведений матриц <math>T_{i}</math>, оно показано на рисунке 1 и использует результаты 1й части.
 +
 
 +
В третьей, последней, части - вычисление результатов, с использованием результатов 2й части.
 +
 
 +
Внутренние графы 1й и 3й частей - пусты, их операции в зависимости только со 2й частью.
 +
 
 +
=== Ресурс параллелизма алгоритма ===
 +
 
 +
На критическом пути алгоритма Стоуна для LU-разложения трёхдиагональной матрицы нужно выполнить:
 +
 
 +
* <math>1</math> ярус делений,
 +
* <math>\lceil \log_2 (n-1) \rceil+1</math> ярусов умножений,
 +
* <math>\lceil \log_2 (n-1) \rceil+1</math> ярусов сложений.
 +
 
 +
Поэтому алгоритм должен быть отнесён к алгоритмам ''логарифмической сложности'' по количеству последовательных операций. Ширина яруса равна <math>n</math>, поэтому алгоритм должен быть отнесён к алгоритмам ''линейной сложности'' по ширине ярусов.
 +
 
 +
=== Входные и выходные данные алгоритма ===
 +
 
 +
'''Входные данные''': трёхдиагональная матрица <math>A</math> (элементы <math>a_{ij}</math>).
 +
 
 +
'''Объём входных данных''': <math>3n-2</math>. В разных реализациях размещение для экономии хранения может быть выполнено разным образом. Например, каждая диагональ может быть строкой массива.
 +
 
 +
'''Выходные данные''': нижняя двухдиагональная матрица <math>L</math> (элементы <math>l_{ij}</math>, причём <math>l_{ii}=1</math>) и верхняя двухдиагональная матрица <math>U</math> (элементы <math>u_{ij}</math>).
 +
 
 +
'''Объём выходных данных''': формально <math>3n-2</math>. Однако благодаря совпадению <math>n</math> входных и выходных данных реально вычисляется лишь <math>2n-2</math> элементов.
 +
 
 +
=== Свойства алгоритма ===
 +
 
 +
Соотношение последовательной и параллельной сложности, как хорошо видно, равно <math>O(n)</math>.
 +
 
 +
При этом вычислительная мощность алгоритма как отношение числа операций к суммарному объему входных и выходных данных является ''логарифмической''.
 +
 
 +
Алгоритм в рамках выбранной версии полностью детерминирован.
 +
 
 +
Устойчивость алгоритма остаётся только при условии невозрастания ведущих главных миноров матрицы. При росте миноров область устойчивости алгоритма очень мала.
 +
 
 +
== Программная реализация алгоритма ==
 +
 
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
−
=== Описание локальности данных и вычислений ===
+
 
−
=== Возможные способы и особенности реализации параллельного алгоритма ===
+
Из-за большой избыточности вычислений алгоритм сдваивания Стоуна никогда не предназначался для последовательной реализации. После обнаружения его неустойчивости стало ясно, что и в будущем он не будет реализован на последовательных архитектурах.
−
=== Масштабируемость алгоритма и его реализации ===
+
 
−
=== Динамические характеристики и эффективность реализации алгоритма ===
+
=== Возможные способы и особенности параллельной реализации алгоритма ===
 +
 
 +
Из-за неустойчивости алгоритма любые способы реализации теряют смысл.
 +
 
 +
Хотя исследование локальности и не имеет смысла из-за неустойчивости, можно сказать, что, как видно по графу алгоритма, ряд дуг длинны как по времени (по различию номеров ярусов операций, являющихся началом и концом дуги), так и по пространству (исключением является только размещение в гиперкубе, физически невозможное). Эта неустранимая нелокальность должна тормозить исполнение алгоритма. Реальное же исследование последовательного кода на обращения в память проводить бессмысленно, поскольку последовательный код не применяется и не будет применяться никем.
 +
 
 +
Из-за неустойчивости алгоритма не имеют смысла любые замеры эффективности любых его реализаций и их динамических характеристик. Кроме того, из-за избыточности вычислений ясно, что реальная эффективность будет существенно ограничена даже при малых размерах задач, где неустойчивость ещё не сказывается.
 +
 
 +
Из-за неустойчивости алгоритма его не используют на практике, поэтому планировавшаяся в исходной публикации<ref name="STONE" /> замена более популярной [[Полный метод циклической редукции|циклической редукции]] не удалась.
 +
 
 +
=== Результаты прогонов ===
 
=== Выводы для классов архитектур ===
 
=== Выводы для классов архитектур ===
−
=== Существующие реализации алгоритма ===
+
 
 +
На области относительно приемлемых ошибок (малые размеры задач) у метода всё равно нет шансов на нормальную реализацию - быстрее будет СЛАУ прогонкой на обычной персоналке решить.
 +
 
 
== Литература ==
 
== Литература ==
  
 
<references />
 
<references />
  
−
[[Категория:Статьи в работе]]
+
[[Категория:Законченные статьи]]
 
[[Категория:Метод сдваивания]]
 
[[Категория:Метод сдваивания]]
 +
[[Категория:Неустойчивые_параллельные_методы]]
 +
[[Категория:Алгоритмы с избыточными вычислениями]]
 +
[[Категория:Разложения матриц]]
 +
 +
[[en:Stone doubling algorithm for the LU decomposition of a tridiagonal matrix]]

Текущая версия на 15:50, 8 июля 2022


Алгоритм сдваивания Стоуна для LU-разложения трёхдиагональной матрицы
Последовательный алгоритм
Последовательная сложность [math]3(n-1)\lceil \log_2 (n-1) \rceil + o((n-1)\lceil \log_2 (n-1) \rceil)[/math]
Объём входных данных [math]3n-2[/math]
Объём выходных данных [math]2n-1[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]2 \lceil \log_2 (n-1) \rceil + 3[/math]
Ширина ярусно-параллельной формы [math]n[/math]


Основные авторы описания: А.В.Фролов

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

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

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

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

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

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

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

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

Для полного выполнения алгоритма Стоуна и получения разложения трёхдиагональной матрицы на две двудиагональные нужно выполнить:

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

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

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

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

В первой части производится вычисление недиагональных элементов матриц [math]T_{i}[/math].

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

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

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

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

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

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

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

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

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

Объём входных данных: [math]3n-2[/math]. В разных реализациях размещение для экономии хранения может быть выполнено разным образом. Например, каждая диагональ может быть строкой массива.

Выходные данные: нижняя двухдиагональная матрица [math]L[/math] (элементы [math]l_{ij}[/math], причём [math]l_{ii}=1[/math]) и верхняя двухдиагональная матрица [math]U[/math] (элементы [math]u_{ij}[/math]).

Объём выходных данных: формально [math]3n-2[/math]. Однако благодаря совпадению [math]n[/math] входных и выходных данных реально вычисляется лишь [math]2n-2[/math] элементов.

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложности, как хорошо видно, равно [math]O(n)[/math].

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

Алгоритм в рамках выбранной версии полностью детерминирован.

Устойчивость алгоритма остаётся только при условии невозрастания ведущих главных миноров матрицы. При росте миноров область устойчивости алгоритма очень мала.

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

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

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

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

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

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

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

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

2.3 Результаты прогонов

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

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

3 Литература

  1. ↑ Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
  2. ↑ Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.
  3. ↑ 3,0 3,1 Stone H.S. An Efficient Parallel Algorithm for the Solution of a Tridiagonal Linear System of Equations // J. ACM, Vol. 20, No. 1 (Jan. 1973), P. 27-38.
Источник — https://algowiki-project.org/w/ru/index.php?title=Алгоритм_сдваивания_Стоуна_для_LU-разложения_трёхдиагональной_матрицы&oldid=29082