Метод сдваивания Стоуна для решения двудиагональных СЛАУ: различия между версиями
Frolov (обсуждение | вклад) |
ASA (обсуждение | вклад) |
||
(не показано 36 промежуточных версий 2 участников) | |||
Строка 1: | Строка 1: | ||
− | {{algorithm}} | + | {{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>, где |
{{Шаблон:Трёхдиагональная СЛАУ в стандартном виде}} | {{Шаблон:Трёхдиагональная СЛАУ в стандартном виде}} | ||
Строка 10: | Строка 19: | ||
[[Метод сдваивания Стоуна]] впервые предложен в начале 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> в качестве альтернативы другим параллельным алгоритмам решения трёхдиагональных СЛАУ, например, [[Полный метод циклической редукции|методу циклической редукции]]. | [[Метод сдваивания Стоуна]] впервые предложен в начале 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> в качестве альтернативы другим параллельным алгоритмам решения трёхдиагональных СЛАУ, например, [[Полный метод циклической редукции|методу циклической редукции]]. | ||
− | Здесь рассматривается его вторая часть - решение двух | + | Здесь рассматривается его вторая часть - решение двух двудиагональных СЛАУ. Оно использует представление матрицы |
{{Шаблон:Трёхдиагональная матрица стандартной формы}} | {{Шаблон:Трёхдиагональная матрица стандартной формы}} | ||
Строка 27: | Строка 36: | ||
=== Математическое описание алгоритма === | === Математическое описание алгоритма === | ||
− | + | Метод Стоуна в части решения двухдиагональных СЛАУ <math>Ly = b</math> и <math>D^{-1}Ux = z</math>, полученных при решении исходной <math>Ax = b</math> после вычисления разложения <math>A = LU</math>, заключается в том, что получающиеся при их непосредственном решении рекурсивные зависимости | |
− | |||
− | [[file:DoublingPartUniver.png|thumb|right|600px|Рисунок | + | [[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_1 = b_1</math>, и потом <math>y_{i} = b_{i} - l_{i i-1} y_{i-1}, i = 2,..., n</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_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>x_{i} = z_{i} - \frac{u_{i i+1}}{u_{ii}} x_{i+1}, i = n-1,...,1</math> | ||
заменяются соответственно на | заменяются соответственно на | ||
Строка 55: | Строка 59: | ||
\end{bmatrix} \begin{bmatrix} | \end{bmatrix} \begin{bmatrix} | ||
y_{i-1} \\ | y_{i-1} \\ | ||
− | 1 \\ \end{bmatrix} = | + | 1 \\ \end{bmatrix} = P_{i} \begin{bmatrix} |
y_{i-1} \\ | y_{i-1} \\ | ||
1 \\ | 1 \\ | ||
− | \end{bmatrix}, i = 2,..., n </math> | + | \end{bmatrix}, i = 2,..., n </math>, где <math> P_{i} = \begin{bmatrix} |
+ | l_{i i-1} & b_{i} \\ | ||
+ | 0 & 1 \\ | ||
+ | \end{bmatrix} | ||
+ | </math> | ||
и | и | ||
Строка 75: | Строка 83: | ||
x_{i+1} \\ | x_{i+1} \\ | ||
1 \\ | 1 \\ | ||
− | \end{bmatrix}, i = n-1,...,1</math> | + | \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> | ||
и после выполнения подстановок оказывается, что | и после выполнения подстановок оказывается, что | ||
Строка 83: | Строка 95: | ||
y_i \\ | y_i \\ | ||
1 \\ | 1 \\ | ||
− | \end{bmatrix} = | + | \end{bmatrix} = P_{i} P_{i-1} ... P_{2} \begin{bmatrix} |
b_{1} \\ | b_{1} \\ | ||
1 \\ | 1 \\ | ||
Строка 99: | Строка 111: | ||
:<math> | :<math> | ||
− | K_{ | + | K_{i} = \begin{bmatrix} |
− | s_{ | + | s_{i} & t_{i} \\ |
0 & 1 \\ | 0 & 1 \\ | ||
− | \end{bmatrix} = | + | \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> | :<math> | ||
Строка 112: | Строка 128: | ||
\end{bmatrix} = C_{i} C_{i+1} ... C_{n-1} , i = 1,...,n-1 </math>. | \end{bmatrix} = C_{i} C_{i+1} ... C_{n-1} , i = 1,...,n-1 </math>. | ||
− | После | + | После чего окончательно вычисляются |
− | |||
− | |||
− | |||
− | |||
<math>x_n = z_n</math>, | <math>x_n = z_n</math>, | ||
Строка 124: | Строка 136: | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
+ | Основную часть вычислительного ядра алгоритма составляют операции трёх типов: деление, умножение и операция типа <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>. | ||
+ | |||
+ | При этом вычислительная мощность алгоритма как отношение числа операций к суммарному объему входных и выходных данных является ''логарифмической''. | ||
+ | |||
+ | Алгоритм в рамках выбранной версии полностью детерминирован. | ||
+ | |||
+ | Для решения СЛАУ с диагональным преобладанием гарантируется устойчивость алгоритма. | ||
== Программная реализация алгоритма == | == Программная реализация алгоритма == | ||
Строка 142: | Строка 206: | ||
Из-за большой избыточности вычислений метод Стоуна никогда не предназначался для последовательной реализации. После обнаружения неустойчивости его первой части стало ясно, что и в будущем он не будет реализован на любых, а не только на последовательных архитектурах. Вторую же часть из-за большой избыточности вычислений тоже нецелесообразно использовать на последовательных компьютерах, где используют, в основном, прогонки разных версий. | Из-за большой избыточности вычислений метод Стоуна никогда не предназначался для последовательной реализации. После обнаружения неустойчивости его первой части стало ясно, что и в будущем он не будет реализован на любых, а не только на последовательных архитектурах. Вторую же часть из-за большой избыточности вычислений тоже нецелесообразно использовать на последовательных компьютерах, где используют, в основном, прогонки разных версий. | ||
− | |||
=== Возможные способы и особенности параллельной реализации алгоритма === | === Возможные способы и особенности параллельной реализации алгоритма === | ||
− | + | ||
− | === | + | В принципе, если взять блочно-двухдиагональный вариант разложения, то, поскольку вес макровершин будет расти больше, чем вес обменов, структуру метода Стоуна вполне можно реализовать с малым временем простоев, но ограничение реальной эффективности обратным двоичным логарифмом делает эти попытки не вполне осмысленными. Видимо, поэтому его и не применяют - ни в точечной формулировке, ни в блочной. |
+ | |||
+ | Как видно по графу алгоритма, ряд дуг длинны как по времени (по различию номеров ярусов операций, являющихся началом и концом дуги), так и по пространству (исключением является только размещение в гиперкубе, физически невозможное). Эта неустранимая нелокальность должна тормозить исполнение алгоритма. Реальное же исследование последовательного кода на обращения в память проводить бессмысленно, поскольку последовательный код не применяется и не будет применяться никем. | ||
+ | |||
+ | При оценке масштабируемости этого алгоритма, как и всех алгоритмов с избыточными вычислениями, следует учитывать, что сравнение по быстродействию и эффективности нужно проводить не с однопроцессорным вариантом исполнения самого алгоритма, а с алгоритмом последовательного решения двух двудиагональных СЛАУ. | ||
+ | Избыточность порядка двоичного логарифма приводит к тому, что реальная эффективность будет уже с размером порядка миллиона примерно 5% от эффективности, достигаемой программистскими ухищрениями, а с дальнейшим ростом процент будет ещё падать. Это соображение ставит крест на попытки эффективно реализовать и масштабировать метод Стоуна. | ||
+ | |||
+ | Из-за большой избыточности вычислений, даже несмотря на устойчивость этой части метода Стоуна, реальная эффективность была бы очень низка. Поэтому его реализации в настоящее время недоступны, и динамические характеристики отсутствуют. | ||
+ | |||
+ | Из-за неустойчивости первой части метода его не используют на практике, поэтому планировавшаяся в исходной публикации<ref name="STONE" /> замена более популярной [[Полный метод циклической редукции|циклической редукции]] не удалась. Реализаций схемы Стоуна отсутствуют в пакетах программ, даже в её второй (устойчивой) части. | ||
+ | |||
+ | === Результаты прогонов === | ||
=== Выводы для классов архитектур === | === Выводы для классов архитектур === | ||
− | |||
− | + | Даже если говорить только о второй части схемы Стоуна, являющейся устойчивой, реальное её использование на любых параллельных архитектурах из-за большой избыточности алгоритма даст очень низкую реальную эффективность. При этом, в отличие от той же циклической редукции, где реальная эффективность того же порядка, из-за большого числа избыточных операций в схеме Стоуна на той же вычислительной системе будет гораздо больше энергопотребление. Поэтому если даже мы имеем полное разложение трёхдиагональной матрицы в произведение двух двудиагональных, вместо схемы Стоуна логичнее для распараллеливания взять алгоритм с существенно большей эффективностью, например последовательно-параллельную схему. | |
== Литература == | == Литература == | ||
Строка 155: | Строка 228: | ||
<references /> | <references /> | ||
− | [[Категория: | + | [[Категория:Законченные статьи]] |
[[Категория:Метод сдваивания]] | [[Категория:Метод сдваивания]] | ||
[[Категория:Алгоритмы с избыточными вычислениями]] | [[Категория:Алгоритмы с избыточными вычислениями]] | ||
− | [[ | + | |
+ | [[en:Stone doubling algorithm for solving bidiagonal SLAEs]] |
Текущая версия на 11:32, 14 июля 2022
Метод сдваивания Стоуна для решения двудиагональных СЛАУ | |
Последовательный алгоритм | |
Последовательная сложность | [math]3(n-1)(\lceil \log_2 (n-1) \rceil+2)[/math] |
Объём входных данных | [math]4n-2[/math] |
Объём выходных данных | [math]n[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math]4 \lceil \log_2 (n-1) \rceil + 5[/math] |
Ширина ярусно-параллельной формы | [math]n[/math] |
Основные авторы описания: А.В.Фролов.
Содержание
- 1 Свойства и структура алгоритмов
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема реализации последовательного алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 3 Литература
1 Свойства и структура алгоритмов
1.1 Общее описание алгоритма
Алгоритм сдваивания Стоуна для решения двудиагональных СЛАУ - часть метода сдваивания Стоуна для решения СЛАУ[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] 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 Математическое описание алгоритма
Метод Стоуна в части решения двухдиагональных СЛАУ [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} = 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].
1.3 Вычислительное ядро алгоритма
Основную часть вычислительного ядра алгоритма составляют операции трёх типов: деление, умножение и операция типа [math]a+bc[/math].
1.4 Макроструктура алгоритма
На макроуровне можно выделить такие 4 макрооперации: вычисление сдваиванием матриц [math]K_{i}[/math], вычисление элементов матриц [math]C_i[/math], вычисление сдваиванием матриц [math]R_{i}[/math], вычисление результатов.
1.5 Схема реализации последовательного алгоритма
Метод Стоуна изначально спроектирован для параллельного исполнения, поскольку является по отношению к, например, классической прогонке, алгоритмом с избыточными вычислениями. Вторая часть метода - для решения двухдиагональных СЛАУ - также содержит избыточные вычисления в сравнении с методами подстановки для их решения.
1.6 Последовательная сложность алгоритма
Для полного выполнения алгоритма Стоуна и решения двух двухдиагональных СЛАУ нужно выполнить (если не используются предвычисления обратных элементов):
- [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] сложений.
Поэтому алгоритм должен быть отнесён к алгоритмам линейно-логарифмической сложности по количеству последовательных операций.
1.7 Информационный граф
Как уже отмечено, макроструктура алгоритма состоит из 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й частей - пусты, их операции в зависимости только с другими частями.
1.8 Ресурс параллелизма алгоритма
На критическом пути алгоритма Стоуна решения двух двухдиагональных СЛАУ нужно выполнить (если не используются предвычисления обратных элементов):
- [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], поэтому алгоритм должен быть отнесён к алгоритмам линейной сложности по ширине ярусов.
1.9 Входные и выходные данные алгоритма
Входные данные: двудиагональные матрицы [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].
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности, как хорошо видно, равно [math]O(n)[/math].
При этом вычислительная мощность алгоритма как отношение числа операций к суммарному объему входных и выходных данных является логарифмической.
Алгоритм в рамках выбранной версии полностью детерминирован.
Для решения СЛАУ с диагональным преобладанием гарантируется устойчивость алгоритма.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
Из-за большой избыточности вычислений метод Стоуна никогда не предназначался для последовательной реализации. После обнаружения неустойчивости его первой части стало ясно, что и в будущем он не будет реализован на любых, а не только на последовательных архитектурах. Вторую же часть из-за большой избыточности вычислений тоже нецелесообразно использовать на последовательных компьютерах, где используют, в основном, прогонки разных версий.
2.2 Возможные способы и особенности параллельной реализации алгоритма
В принципе, если взять блочно-двухдиагональный вариант разложения, то, поскольку вес макровершин будет расти больше, чем вес обменов, структуру метода Стоуна вполне можно реализовать с малым временем простоев, но ограничение реальной эффективности обратным двоичным логарифмом делает эти попытки не вполне осмысленными. Видимо, поэтому его и не применяют - ни в точечной формулировке, ни в блочной.
Как видно по графу алгоритма, ряд дуг длинны как по времени (по различию номеров ярусов операций, являющихся началом и концом дуги), так и по пространству (исключением является только размещение в гиперкубе, физически невозможное). Эта неустранимая нелокальность должна тормозить исполнение алгоритма. Реальное же исследование последовательного кода на обращения в память проводить бессмысленно, поскольку последовательный код не применяется и не будет применяться никем.
При оценке масштабируемости этого алгоритма, как и всех алгоритмов с избыточными вычислениями, следует учитывать, что сравнение по быстродействию и эффективности нужно проводить не с однопроцессорным вариантом исполнения самого алгоритма, а с алгоритмом последовательного решения двух двудиагональных СЛАУ. Избыточность порядка двоичного логарифма приводит к тому, что реальная эффективность будет уже с размером порядка миллиона примерно 5% от эффективности, достигаемой программистскими ухищрениями, а с дальнейшим ростом процент будет ещё падать. Это соображение ставит крест на попытки эффективно реализовать и масштабировать метод Стоуна.
Из-за большой избыточности вычислений, даже несмотря на устойчивость этой части метода Стоуна, реальная эффективность была бы очень низка. Поэтому его реализации в настоящее время недоступны, и динамические характеристики отсутствуют.
Из-за неустойчивости первой части метода его не используют на практике, поэтому планировавшаяся в исходной публикации[3] замена более популярной циклической редукции не удалась. Реализаций схемы Стоуна отсутствуют в пакетах программ, даже в её второй (устойчивой) части.
2.3 Результаты прогонов
2.4 Выводы для классов архитектур
Даже если говорить только о второй части схемы Стоуна, являющейся устойчивой, реальное её использование на любых параллельных архитектурах из-за большой избыточности алгоритма даст очень низкую реальную эффективность. При этом, в отличие от той же циклической редукции, где реальная эффективность того же порядка, из-за большого числа избыточных операций в схеме Стоуна на той же вычислительной системе будет гораздо больше энергопотребление. Поэтому если даже мы имеем полное разложение трёхдиагональной матрицы в произведение двух двудиагональных, вместо схемы Стоуна логичнее для распараллеливания взять алгоритм с существенно большей эффективностью, например последовательно-параллельную схему.
3 Литература
- ↑ Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
- ↑ Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.
- ↑ 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.