Метод сдваивания Стоуна для решения двудиагональных СЛАУ
| Метод сдваивания Стоуна для решения двудиагональных СЛАУ | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]\displaystyle{ 3(n-1)(\lceil \log_2 (n-1) \rceil+2) }[/math] |
| Объём входных данных | [math]\displaystyle{ 4n-2 }[/math] |
| Объём выходных данных | [math]\displaystyle{ n }[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]\displaystyle{ 4 \lceil \log_2 (n-1) \rceil + 5 }[/math] |
| Ширина ярусно-параллельной формы | [math]\displaystyle{ n }[/math] |
Основные авторы описания: А.В.Фролов
Свойства и структура алгоритмов
Общее описание алгоритма
Алгоритм сдваивания Стоуна для решения двудиагональных СЛАУ - часть метода сдваивания Стоуна для решения СЛАУ[1][2] вида [math]\displaystyle{ Ax = b }[/math], где
- [math]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ 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]\displaystyle{ Ax = b }[/math] можно поменять на последовательное решение двух СЛАУ [math]\displaystyle{ Ly = b }[/math] и затем [math]\displaystyle{ Ux = y }[/math]. При этом вторую СЛАУ тоже можно решить как последовательность СЛАУ [math]\displaystyle{ Dz = y }[/math] и [math]\displaystyle{ D^{-1}Ux = z }[/math], где [math]\displaystyle{ D }[/math] - диагональная матрица, составленная из диагональных элементов матрицы [math]\displaystyle{ U }[/math].
Математическое описание алгоритма
Метод Стоуна в части решения двухдиагональных СЛАУ [math]\displaystyle{ Ly = b }[/math] и [math]\displaystyle{ D^{-1}Ux = z }[/math], полученных при решении исходной [math]\displaystyle{ Ax = b }[/math] после вычисления разложения [math]\displaystyle{ A = LU }[/math], заключается в том, что получающиеся при их непосредственном решении рекурсивные зависимости
[math]\displaystyle{ y_1 = b_1 }[/math], и потом [math]\displaystyle{ y_{i} = b_{i} - l_{i i-1} y_{i-1}, i = 2,..., n }[/math]
и
[math]\displaystyle{ x_n = z_n }[/math], и потом [math]\displaystyle{ x_{i} = z_{i} - \frac{u_{i i+1}}{u_{ii}} x_{i+1}, i = n-1,...,1 }[/math]
заменяются соответственно на
- [math]\displaystyle{ \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]\displaystyle{ P_{i} = \begin{bmatrix} l_{i i-1} & b_{i} \\ 0 & 1 \\ \end{bmatrix} }[/math]
и
- [math]\displaystyle{ \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]\displaystyle{ C_{i} = \begin{bmatrix} \frac{u_{i i+1}}{u_{ii}} & z_{i} \\ 0 & 1 \\ \end{bmatrix} }[/math]
и после выполнения подстановок оказывается, что
- [math]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ K_{i} = \begin{bmatrix} s_{i} & t_{i} \\ 0 & 1 \\ \end{bmatrix} = P_{i} P_{n-1} ... P_{2} , i = 2,..., n }[/math]
После этого вычисляются промежуточные результаты
[math]\displaystyle{ z_i = \frac{s_{i}b_{1}+t_{i}}{u_{ii}}, i = 2,..., n }[/math],
и снова по схеме сдваивания матрицы
- [math]\displaystyle{ 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]\displaystyle{ x_n = z_n }[/math],
[math]\displaystyle{ x_i = v_{i}z_{n}+w_{i}, i = 1,...,n-1 }[/math].
Вычислительное ядро алгоритма
Основную часть вычислительного ядра алгоритма составляют операции трёх типов: деление, умножение и операция типа [math]\displaystyle{ a+bc }[/math].
Макроструктура алгоритма
На макроуровне можно выделить такие 4 макрооперации: вычисление сдваиванием матриц [math]\displaystyle{ K_{i} }[/math], вычисление элементов матриц [math]\displaystyle{ C_i }[/math], вычисление сдваиванием матриц [math]\displaystyle{ R_{i} }[/math], вычисление результатов.
Схема реализации последовательного алгоритма
Метод Стоуна изначально спроектирован для параллельного исполнения, поскольку является по отношению к, например, классической прогонке, алгоритмом с избыточными вычислениями. Вторая часть метода - для решения двухдиагональных СЛАУ - также содержит избыточные вычисления в сравнении с методами подстановки для их решения.
Последовательная сложность алгоритма
Для полного выполнения алгоритма Стоуна и решения двух двухдиагональных СЛАУ нужно выполнить (если не используются предвычисления обратных элементов):
- [math]\displaystyle{ 2n-2 }[/math] делений,
- [math]\displaystyle{ 2(n-1)(\lceil \log_2 (n-1) \rceil+1) }[/math] умножений,
- [math]\displaystyle{ (n-1)(\lceil \log_2 (n-1) \rceil+2) }[/math] сложений.
Поэтому алгоритм должен быть отнесён к алгоритмам линейно-логарифмической сложности по количеству последовательных операций.
Информационный граф
Как уже отмечено, макроструктура алгоритма состоит из 4 частей.
В первой части производится вычисление сдваиванием матриц [math]\displaystyle{ K_{i} }[/math]. Граф этой части показан на Рисунке 1.
Во второй части происходит вычисление элементов матриц [math]\displaystyle{ C_i }[/math]. При этом внедиагональные элементы ([math]\displaystyle{ z_i }[/math]) вычисляются с использованием результатов первой части, а вот диагональные (частные [math]\displaystyle{ \frac{u_{i i+1}}{u_{ii}} }[/math]) можно вычислить заранее, одновременно с 1й частью.
В третьей части - вычисление сдваиванием матриц [math]\displaystyle{ R_{i} }[/math], оно показано на рисунке 2 и использует результаты 2й части.
В четвёртой, последней, части - вычисление результатов, с использованием результатов 2й части.
Внутренние графы 2й и 4й частей - пусты, их операции в зависимости только с другими частями.
Ресурс параллелизма алгоритма
На критическом пути алгоритма Стоуна решения двух двухдиагональных СЛАУ нужно выполнить (если не используются предвычисления обратных элементов):
- [math]\displaystyle{ 1 }[/math] ярус делений,
- [math]\displaystyle{ 2(\lceil \log_2 (n-1) \rceil+2) }[/math] ярусов умножений,
- [math]\displaystyle{ 2(\lceil \log_2 (n-1) \rceil+2) }[/math] ярусов сложений.
Поэтому алгоритм должен быть отнесён к алгоритмам логарифмической сложности по количеству последовательных операций. Ширина яруса равна [math]\displaystyle{ n }[/math], поэтому алгоритм должен быть отнесён к алгоритмам линейной сложности по ширине ярусов.
Входные и выходные данные алгоритма
Входные данные: двудиагональные матрицы [math]\displaystyle{ L }[/math] (с единичной диагональю) (элементы [math]\displaystyle{ l_{ij} }[/math]) и [math]\displaystyle{ U }[/math] (элементы [math]\displaystyle{ u_{ij} }[/math]), вектор [math]\displaystyle{ b }[/math] (элементы [math]\displaystyle{ b_{i} }[/math]).
Объём входных данных: [math]\displaystyle{ 4n-2 }[/math].
Выходные данные: вектор [math]\displaystyle{ x }[/math] (элементы [math]\displaystyle{ x_{i} }[/math]).
Объём выходных данных: [math]\displaystyle{ n }[/math].
Свойства алгоритма
Соотношение последовательной и параллельной сложности, как хорошо видно, равно [math]\displaystyle{ O(n) }[/math].
При этом вычислительная мощность алгоритма как отношение числа операций к суммарному объему входных и выходных данных является логарифмической.
Алгоритм в рамках выбранной версии полностью детерминирован.
Для решения СЛАУ с диагональным преобладанием гарантируется устойчивость алгоритма.
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Из-за большой избыточности вычислений метод Стоуна никогда не предназначался для последовательной реализации. После обнаружения неустойчивости его первой части стало ясно, что и в будущем он не будет реализован на любых, а не только на последовательных архитектурах. Вторую же часть из-за большой избыточности вычислений тоже нецелесообразно использовать на последовательных компьютерах, где используют, в основном, прогонки разных версий.
Локальность данных и вычислений
Как видно по графу алгоритма, ряд дуг длинны как по времени (по различию номеров ярусов операций, являющихся началом и концом дуги), так и по пространству (исключением является только размещение в гиперкубе, физически невозможное). Эта неустранимая нелокальность должна тормозить исполнение алгоритма. Реальное же исследование последовательного кода на обращения в память проводить бессмысленно, поскольку последовательный код не применяется и не будет применяться никем.
Возможные способы и особенности параллельной реализации алгоритма
В принципе, если взять блочно-двухдиагональный вариант разложения, то, поскольку вес макровершин будет расти больше, чем вес обменов, структуру метода Стоуна вполне можно реализовать с малым временем простоев, но ограничение реальной эффективности обратным двоичным логарифмом делает эти попытки не вполне осмысленными. Видимо, поэтому его и не применяют - ни в точечной формулировке, ни в блочной.
Масштабируемость алгоритма и его реализации
При оценке масштабируемости этого алгоритма, как и всех алгоритмов с избыточными вычислениями, следует учитывать, что сравнение по быстродействию и эффективности нужно проводить не с однопроцессорным вариантом исполнения самого алгоритма, а с алгоритмом последовательного решения двух двудиагональных СЛАУ. Избыточность порядка двоичного логарифма приводит к тому, что реальная эффективность будет уже с размером порядка миллиона примерно 5% от эффективности, достигаемой программистскими ухищрениями, а с дальнейшим ростом процент будет ещё падать. Это соображение ставит крест на попытки эффективно реализовать и масштабировать метод Стоуна.
Динамические характеристики и эффективность реализации алгоритма
Из-за большой избыточности вычислений, даже несмотря на устойчивость этой части метода Стоуна, реальная эффективность была бы очень низка. Поэтому его реализации в настоящее время недоступны, и динамические характеристики отсутствуют.
Выводы для классов архитектур
Даже если говорить только о второй части схемы Стоуна, являющейся устойчивой, реальное её использование на любых параллельных архитектурах из-за большой избыточности алгоритма даст очень низкую реальную эффективность. При этом, в отличие от той же циклической редукции, где реальная эффективность того же порядка, из-за большого числа избыточных операций в схеме Стоуна на той же вычислительной системе будет гораздо больше энергопотребление. Поэтому если даже мы имеем полное разложение трёхдиагональной матрицы в произведение двух двудиагональных, вместо схемы Стоуна логичнее для распараллеливания взять алгоритм с существенно большей эффективностью, например последовательно-параллельную схему.
Существующие реализации алгоритма
Из-за неустойчивости первой части метода его не используют на практике, поэтому планировавшаяся в исходной публикации[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.