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

{{algorithm
| name              = Метод сдваивания Стоуна для решения двудиагональных СЛАУ
| serial_complexity = <math>3(n-1)(\lceil \log_2 (n-1) \rceil+2)</math>
| pf_height         = <math>4 \lceil \log_2 (n-1) \rceil + 5</math>
| pf_width          = <math>n</math>
| input_data        = <math>4n-2</math>
| output_data       = <math>n</math>
}}

Основные авторы описания: [[Участник:Frolov|А.В.Фролов]].

== Свойства и структура алгоритмов ==
=== Общее описание алгоритма ===

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

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

{{Шаблон:Трёхдиагональная матрица стандартной формы}}

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

{{Шаблон:L2dUniDiag}}

и 

{{Шаблон:U2dCommon}}

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

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

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

[[file:DoublingPartUniver.png|thumb|right|600px|Рисунок 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>.

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

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

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

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

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

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

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

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

* <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> сложений. 

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

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

Как уже отмечено, макроструктура алгоритма состоит из 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й частей - пусты, их операции в зависимости только с другими частями.

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

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

* <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>, поэтому алгоритм должен быть отнесён к алгоритмам ''линейной сложности'' по ширине ярусов.

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

'''Входные данные''': двудиагональные матрицы <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>.

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

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

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

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

Для решения СЛАУ с диагональным преобладанием гарантируется устойчивость алгоритма.

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

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

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

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

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

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

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

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

Из-за неустойчивости первой части метода его не используют на практике, поэтому планировавшаяся в исходной публикации<ref name="STONE" /> замена более популярной [[Полный метод циклической редукции|циклической редукции]] не удалась. Реализаций схемы Стоуна отсутствуют в пакетах программ, даже в её второй (устойчивой) части. 

=== Результаты прогонов ===
=== Выводы для классов архитектур ===

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

== Литература ==

<references />

[[Категория:Законченные статьи]]
[[Категория:Метод сдваивания]]
[[Категория:Алгоритмы с избыточными вычислениями]]

[[en:Stone doubling algorithm for solving bidiagonal SLAEs]]