Elimination method, pointwise version

Primary authors of this description: A.V.Frolov, Vad.V.Voevodin (Section 2.2), A.M.Teplov (Section 2.4)

Contents

1 Properties and structure of the algorithm 1.1 General description of the algorithm 1.2 Mathematical description of the algorithm 1.3 Computational kernel of the algorithm 1.4 Macro structure of the algorithm 1.5 Implementation scheme of the serial algorithm 1.6 Serial complexity of the algorithm 1.7 Information graph 1.8 Parallelization resource of the algorithm 1.9 Input and output data of the algorithm 1.10 Properties of the algorithm 2 Software implementation of the algorithm 2.1 Implementation peculiarities of the serial algorithm 2.2 Locality of data and computations 2.2.1 Locality of implementation 2.2.1.1 Structure of memory access and a qualitative estimation of locality 2.2.1.2 Quantitative estimation of locality 2.3 Possible methods and considerations for parallel implementation of the algorithm 2.4 Scalability of the algorithm and its implementations 2.4.1 Scalability of the algorithm 2.4.2 Scalability of of the algorithm implementation 2.5 Dynamic characteristics and efficiency of the algorithm implementation 2.6 Conclusions for different classes of computer architecture 2.7 Existing implementations of the algorithm 3 References

1 Properties and structure of the algorithm

1.1 General description of the algorithm

The elimination method is a variant of Gaussian elimination used for solving a system of linear algebraic equations SLAE) [1][2] Ax = b, where

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}

However, presentations of the elimination method [3] often use a different notation and a numbering for the right-hand side and matrix of the system. For instance, the above SLAE can be written as

A = \begin{bmatrix}

c_{0} & -b_{0}  & 0      &  \cdots & \cdots & 0 \\
-a_{1} & c_{1}  & -b_{1} &  \cdots & \cdots & 0 \\
0      &  -a_{2} & c_{2} &  \cdots & \cdots & 0 \\
\vdots & \vdots & \ddots & \ddots & \ddots & 0 \\
0 & \cdots & \cdots & -a_{N-1} & c_{N-1}  & -b_{N-1} \\
0 & \cdots & \cdots & 0 & -a_{N}  & c_{N} \\

\end{bmatrix}\begin{bmatrix} y_{0} \\ y_{1} \\ \vdots \\ y_{N} \\ \end{bmatrix} = \begin{bmatrix} f_{0} \\ f_{1} \\ \vdots \\ f_{N} \\ \end{bmatrix}

(here, N=n-1). If each equation is written separately, then we have

c_{0} y_{0} - b_{0} y_{1} = f_{0}, -a_{i} y_{i-1} + c_{i} y_{i} - b_{i} y_{i+1} = f_{i}, 1 \le i \le N-1, -a_{N} y_{N-1} + c_{N} y_{N} = f_{N}.

Here, we examine the so-called right-elimination, which is the variant of elimination method in which a SLAE is processed top-down and then in the reverse direction. In this variant, one first goes top down eliminating the subdiagonal unknowns and then bottom up eliminating the superdiagonal ones. There is also a variant, called the left-elimination, in which a SLAE is first processed bottom-up and then top-down. There is basically no difference between both variants; consequently, we do not include a separate description of the left-elimination.

Figure 1. Graph of the elimination method for n=8. The input and output data are not shown. The symbols / and f denote division and operation a+bc or a-bc, respectively.

1.2 Mathematical description of the algorithm

In the notation introduced above, the forward elimination path consists in calculating the elimination coefficients

\alpha_{1} = b_{0}/c_{0}, \beta_{1} = f_{0}/c_{0}, \alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i}), \quad i = 1, 2, \cdots , N-1, \beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i}), \quad i = 1, 2, \cdots , N.

Then, in the backward elimination path, one calculates the solution

y_{N} = \beta_{N+1}, y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}, \quad i = N-1, N-2, \cdots , 1, 0.

According to the literature (see [3]), these formulas are equivalent to calculating a variant of the LU decomposition for the coefficient matrix accompanied by solving bidiagonal systems via forward and back substitutions.

1.3 Computational kernel of the algorithm

The computational kernel of this algorithm can be thought of as compiled of two parts, namely, the forward and backward elimination paths. The computational kernel of the forward elimination path consists of sequences of divisions, multiplications, and additions/subtractions. The computational kernel of the backward elimination path contains only multiplications and additions sequences.

Figure 2. Detailed graph of the elimination method for n=8 (each reciprocal number is only once calculated). The input and output data are not shown. Each calculation of a reciprocal number is marked by inv, and each multiplication is marked by mult. The operations repeated after the replacement of the right-hand side of SLAE are set out in grey.

1.4 Macro structure of the algorithm

The macro structure of this algorithm can be represented as the combination of the forward and backward elimination paths. In addition, the forward path can also be split into two macro units, namely, the triangular decomposition of the coefficient matrix and the forward substitution for a bidiagonal SLAE. These units are executed simultaneously, that is, in parallel. The process of solving the bidiagonal SLAE uses the results produced by the triangular decomposition.

1.5 Implementation scheme of the serial algorithm

The method is executed as the following sequence of steps:

1. Initialize the forward elimination path:

\alpha_{1} = b_{0}/c_{0}, \beta_{1} = f_{0}/c_{0}

2. For i increasing from 1 to N-1, execute the forward path formulas:

\alpha_{i+1} = b_{i}/(c_{i}-a_{i}\alpha_{i}), \beta_{i+1} = (f_{i}+a_{i}\beta_{i})/(c_{i}-a_{i}\alpha_{i}).

3. Initialize the backward elimination path:

y_{N} = (f_{N}+a_{N}\beta_{N})/(c_{N}-a_{N}\alpha_{N})

4. For i decreasing from N-1 to 0, execute the backward path formulas: y_{i} = \alpha_{i+1} y_{i+1} + \beta_{i+1}.

The formulas of the forward path contain double divisions by the same expressions. These divisions can be replaced by the calculation of reciprocal numbers succeeded by the multiplication by these numbers.

1.6 Serial complexity of the algorithm

Consider a tridiagonal SLAE consisting of n equations with n unknowns. The elimination method as applied to solving such a SLAE in its (fastest) serial form requires 2n-1 divisions, 3n-3 additions/subtractions, and 3n-3 multiplications. Thus, in terms of serial complexity, the elimination method is qualified as a linear complexity algorithm.

1.7 Information graph

The information graph of the elimination method is shown in fig. 1. An analysis indicates that this graph is virtually serial. During the forward path, two branches (namely, the left branch, corresponding to the matrix decomposition, and the central branch, responsible for the solution of the first bidiagonal system) can be executed in parallel. The right branch corresponds to the backward path. It is evident from the figure that not only the mathematical content of this process but also the structure of the algorithm graph and the direction of the corresponding data flows are in complete agreement with the appellation "backward path." The variant with divisions replaced by reciprocal number calculations is illustrated by the graph in fig. 2.

1.8 Parallelization resource of the algorithm

The parallel version of the elimination method as applied to solving a tridiagonal SLAE consisting of n equations with n unknowns requires that the following tiers be executed:

n division tiers (all the tiers except for the one contain 2 divisions),

2n - 2 multiplication tiers and the same number of addition/subtraction tiers (n-1 tiers contain two operations of both types, while another n-1 tiers contain a single operation of each type).

Thus, in terms of the parallel form height, the elimination method is qualified as an algorithm of complexity O(n). In terms of the parallel form width, its complexity is 2.

1.9 Input and output data of the algorithm

Input data: tridiagonal matrix A (with entries a_{ij}), vector b (with components b_{i}).

Output data: vector x (with components x_{i}).

Size of the output data: n.

1.10 Свойства алгоритма[править] Соотношение последовательной и параллельной сложности, как хорошо видно, является константой (причём менее 2). При этом вычислительная мощность алгоритма как отношение числа операций к суммарному объему входных и выходных данных также является константой. Алгоритм в рамках выбранной версии полностью детерминирован. Обычно прогонка используется для решения СЛАУ с диагональным преобладанием. В этом случае гарантируется устойчивость алгоритма. В случае, когда требуется решение нескольких СЛАУ с одной и той же матрицей, левую ветвь вычислений (см. рисунки с графом алгоритма) можно не повторять. Это связано с тем, что LU-разложение матрицы системы не нужно перевычислять. Тогда будет более предпочтителен вариант с заменой делений. 2 Программная реализация алгоритма[править] 2.1 Особенности реализации последовательного алгоритма[править] В зависимости от способа хранения матрицы СЛАУ (в виде одного массива с 3 строками или в виде 3 разных массивов) и способа хранения вычисляемых коэффициентов (на месте уже использованных элементов матрицы либо отдельно) возможны различные реализации алгоритма. Приведем пример подпрограммы на языке Fortran, реализующей прогонку, где все элементы матрицы хранятся в одном массиве, причём соседние элементы матричной строки размещаются рядом, а вычисляемые коэффициенты --- на месте элементов исходной матрицы.

     subroutine progm (a,x,N)
     real a(3,0:N), x(0:N)
     a(2,0)=1./a(2,0)
     a(3,0)=-a(3,0)*a(2,0) ! alpha 1
     x(0)=x(0)*a(2,0) ! beta 1
     do 10 i=1,N-1
       a(2,i)=1./(a(2,i)+a(1,i)*a(2,i-1))
       a(3,i)=-a(3,i)*a(2,i) ! alpha i+1
       x(i)=(x(i)-a(1,i)*x(i-1))*a(2,i) ! beta i+1
 10  continue
     a(2,N)=1./(a(2,N)+a(1,N)*a(2,N-1))
     x(N)=(x(N)-a(1,N)*x(N-1))*a(2,N) ! y N
     do 20 i=N-1,0,-1
       x(i)=a(3,i)*x(i+1)+x(i) ! y i
 20  continue
     return
     end

В Сети имеется много простых реализаций метода, например в Викиучебник - Реализации алгоритмов - Метод прогонки 2.2 Локальность данных и вычислений[править] Как видно по графу алгоритма, локальность данных по пространству хорошая - все аргументы, которые необходимы для операций, вычисляются "рядом". Однако по времени локальность вычислений не столь хороша. Если данные задачи не помещаются в кэш, то вычисления в "верхнем левом углу" СЛАУ будут выполняться с постоянными промахами кэша. Отсюда может следовать одна из рекомендаций прикладникам, использующим прогонку, - нужно организовать все вычисления так, чтобы прогонки были "достаточно коротки" для помещения данных в кэш. При этом, однако, прогонка относится к таким последовательным алгоритмам, в которых локальность вычислений настолько велика, что является даже излишней[4]. Из-за того, что данные, необходимые для выполнения основных операций прогонки, вычисляются в непосредственно предшествующим им операциях, возможность использования суперскалярности вычислительных ядер процессоров практически сводится на нет, что резко ухудшает эффективность выполнения прогонки даже на современных однопроцессорных и одноядерных системах. 2.2.1 Локальность реализации алгоритма[править] 2.2.1.1 Структура обращений в память и качественная оценка локальности[править]

Рисунок 3. Прогонка, точечный вариант. Общий профиль обращений в память На рис. 3 представлен профиль обращений в память для реализации точечного варианта прогонки. Данный профиль состоит из обращений к 5 массивам – трем диагоналям матрицы, вектору правых частей и результирующему вектору. Исходя из общего профиля, можно увидеть, что программа состоит из двух этапов – на первом выполняется последовательный перебор элементов 4-х массивов, на втором – также последовательный перебор, только в обратном порядке.

Рисунок 4. Фрагмент общего профиля (выделенный зеленым цветом на рис. 3) Однако при более детальном рассмотрении можно заметить, что профиль устроен несколько сложнее (рис. 4). Для некоторых массивов на каждом шаге выполняется по несколько очень близких обращений, что приводит к улучшению как пространственной, так и временной локальности. Таким образом, общий профиль обладает высокой пространственной локальностью (поскольку здесь доминируют последовательные переборы элементов) и средней временной (поскольку некоторые элементы сразу используются повторно). 2.2.1.2 Количественная оценка локальности[править] Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен здесь (функция Kernel). Условия запуска описаны здесь. Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg. На рисунке 5 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Для данной программы значения daps является достаточно неожиданным. Согласно качественному анализу, локальность обращений в память в этом случае достаточно высока. Однако daps оценивает производительность работы с памятью ниже, чем для теста RandomAccess! Такое различие между локальностью и производительностью встречается нечасто. При детальном рассмотрении исходного кода выяснились две причины такого поведения. Во-первых, итерации основного цикла являются зависимыми по памяти, то есть на очередной итерации используются данные, посчитанные на предыдущей. Похоже, что в данном случае эта зависимость значительно снижает эффективность работы с памятью, в то время как с точки зрения локальности данный случай является достаточно неплохим. Также выяснилось, что одна из операций деления неожиданно существенно замедляет выполнение программы. При замене этой операции на умножение время выполнения уменьшается в 2.5 раза, при этом ассемблерный код двух вариантов программы выглядел идентично (за исключением очевидной замены одной команды деления на команду умножения).

Рисунок 5. Сравнение значений оценки daps Такие особенности программы не отражаются в ее свойствах локальности, что и привело к такому существенному разрыву между оценкой локальности и производительностью работы с памятью. Вторая характеристика – cvg – предназначена для получения более машинно-независимой оценки локальности. Она определяет, насколько часто в программе необходимо подтягивать данные в кэш-память. Соответственно, чем меньше значение cvg, тем реже это нужно делать, и тем лучше локальность алгоритма. На рисунке 6 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). В данном случае все предыдущие выводы подтверждаются – локальность оказывается достаточно высокой, что и было замечено согласно качественной оценке локальности.

Рисунок 6. Сравнение значений оценки cvg 2.3 Возможные способы и особенности параллельной реализации алгоритма[править] Как видно из анализа графа алгоритма, его (без существенных модификаций) практически невозможно распараллелить. Поэтому есть два способа использования прогонок для параллельных вычислительных систем: либо разбивать задачу, где используются прогонки, так, чтобы их было достаточно много, например, так, чтобы на каждую из прогонок приходился 1 процессор (1 ядро), либо использовать вместо прогонки её параллельные варианты (циклическую редукцию, последовательно-параллельные варианты и т.п.). 2.4 Масштабируемость алгоритма и его реализации[править] О масштабируемости самой прогонки, как полностью непараллельного алгоритма, говорить не имеет смысла. Однако необходимо отметить, что анализ масштабируемости параллельных вариантов прогонки должен проводиться относительно однопроцессорной реализации описанного классического варианта прогонки, а не относительно однопроцессорных расчетов для её параллельных вариантов. Проведём исследование масштабируемости реализации алгоритма согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов"[5] Суперкомпьютерного комплекса Московского университета. Набор и границы значений изменяемых параметров запуска реализации алгоритма: число процессоров 1; размер области [10240 : 1024000] с шагом 10240. В результате проведенных экспериментов был получен следующий диапазон эффективности реализации алгоритма: минимальная эффективность реализации 0.019%; максимальная эффективность реализации 0.0255%. На следующих рисунках приведены графики производительности и эффективности выбранной реализации алгоритма в зависимости от изменяемых параметров запуска.

Рисунок 7. Реализация алгоритма. Изменение производительности в зависимости от размера вектора.

Рисунок 8. Реализация алгоритма. Изменение эффективности в зависимости от размера вектора. Малая эффективность, по-видимому, связана с избыточной локальностью, описанной в разделе о локальности данных и вычислений. Исследованная реализация алгоритма 2.5 Динамические характеристики и эффективность реализации алгоритма[править] В силу существенно последовательной природы алгоритма и его избыточной локальности, исследование его динамических характеристик представляется малоценным и потому не проводилось. 2.6 Выводы для классов архитектур[править] Прогонка - метод для архитектуры классического, фон-неймановского типа, непригодный даже для эффективной загрузки одноядерных систем, поддерживающих суперскалярность. Для распараллеливания решения СЛАУ с трёхдиагональной матрицей следует использовать какой-либо её параллельный заменитель, например, наиболее распространённую циклическую редукцию, или уступающий ей по критическому пути графа, но имеющий более регулярную структуру графа новый последовательно-параллельный вариант. 2.7 Существующие реализации алгоритма[править] Алгоритм прогонки является настолько простым, что, несмотря на его "дежурное" присутствие в стандартных пакетах программ решения задач линейной алгебры, большинство использующих его исследователей-прикладников часто пишут соответствующий фрагмент программы самостоятельно. Однако надо отметить, что, как правило, метод прогонки в пакетах реализуют не в его "чистом виде", а в виде пары "разложение на две двухдиагональные матрицы" и "решение СЛАУ с предварительно факторизованной трёхдиагональной матрицей". Так обстоит дело, например, в пакете SCALAPACK и его предшественниках. 3 Литература[править] Перейти ↑ Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977. Перейти ↑ Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984. ↑ Перейти к: 3,0 3,1 Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. Перейти ↑ Фролов А.В., Антонов А.С., Воеводин Вл.В., Теплов А.М. Сопоставление разных методов решения одной задачи по методике проекта Algowiki // Подана на конференцию "Параллельные вычислительные технологии (ПаВТ'2016)". Перейти ↑ Воеводин Вл., Жуматий С., Соболев С., Антонов А., Брызгалов П., Никитенко Д., Стефанов К., Воеводин Вад. Практика суперкомпьютера «Ломоносов» // Открытые системы, 2012, N 7, С. 36-39. Категории: Уровень алгоритмаЗаконченные статьиАлгоритмы с низким уровнем параллелизма Навигация ИкрамовОбсуждениеНастройкиСписок наблюденияВкладВыйтиСтатьяОбсуждениеЧитатьНепроверенные измененияПравитьИсторияНе следить

Поиск Перейти Заглавная страница Общий форум Технический форум Справка Свежие правки Хранилище файлов Новые файлы Загрузить файл Инструменты Ссылки сюда Связанные правки Загрузить файл Спецстраницы Версия для печати Постоянная ссылка Сведения о странице На других языках English Последнее изменение этой страницы: 10:08, 28 января 2016. Содержимое доступно по лицензии Creative Commons Attribution