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

Полный метод циклической редукции

Материал из Алговики
Перейти к навигации Перейти к поиску


Циклическая редукция для трёхдиагональной матрицы,
точечный вариант
Последовательный алгоритм
Последовательная сложность [math]17n + o(n)[/math]
Объём входных данных [math]4n-2[/math]
Объём выходных данных [math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]7 log_2 n[/math]
Ширина ярусно-параллельной формы [math]3n/2[/math]


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

Содержание

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]

Бывает, однако, что при изложении сути методов решения трёхдиагональных СЛАУ[3] элементы правой части и матрицы системы обозначают и нумеруют по-другому, например СЛАУ может иметь вид

[math] A = \begin{bmatrix} b_{1} & a_{1} & 0 & \cdots & \cdots & 0 \\ c_{2} & b_{2} & a_{2} & \cdots & \cdots & 0 \\ 0 & c_{3} & b_{3} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & c_{n-1} & b_{n-1} & a_{n-1} \\ 0 & \cdots & \cdots & 0 & c_{n} & b_{n} \\ \end{bmatrix}\begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \\ \end{bmatrix} = \begin{bmatrix} f_{1} \\ f_{2} \\ \vdots \\ f_{n} \\ \end{bmatrix} [/math]

или, если записывать отдельно по уравнениям, то

[math]b_{1} x_{1} + a_{1} x_{2} = f_{1}[/math],

[math]c_{i} x_{i-1} + b_{i} x_{i} + a_{i} x_{i+1} = f_{i}, 2 \le i \le n-1[/math],

[math]c_{n} x_{n-1} + b_{n} x_{n} = f_{n}[/math].

Циклическая редукция, как и все варианты прогонки, заключается [3][4] в исключении из уравнений неизвестных, однако, в отличие от них, в ней исключение ведут одновременно по всей СЛАУ. В принципе, её можно считать вариантом метода редукции, выполняемого максимально возможное для данной СЛАУ число раз.

1.2 Математическое описание алгоритма

Лучше всего схема циклической редукции[3] разработана для случая [math]n = 2^{m}-1[/math]. Эта схема состоит из прямого и обратного ходов. Прямой ход состоит из последовательного уменьшения в СЛАУ количества уравнений почти в 2 раза (за счёт подстановки из уравнений с нечётными номерами заменяются уравнения с чётными), пока не останется одно уравнение, обратный - в получении всё большего количества компонент решения исходной СЛАУ. Оба хода - как прямой, так и обратный - разбиты на шаги. Здесь мы приведём тот вариант алгоритма, в котором операции экономятся за счёт предварительной нормировки уравнений, используемых для исключения неизвестных.

1.2.1 Прямой ход редукции

В начале считается, что все [math]c^{(0)}_{i} = c_{i}, b^{(0)}_{i} = b_{i}, a^{(0)}_{i} = a_{i}, f^{(0)}_{i} = f_{i}, x^{(0)}_{i} = x_{i}[/math]

На k-м шаге теперь выполняем процедуру редукции системы уравнений размерности n.

Для каждого из уравнений

[math]c^{(k)}_{i} x^{(k)}_{i-1} + b^{(k)}_{i} x^{(k)}_{i} + a^{(k)}_{i} x^{(k)}_{i+1} = f^{(k)}_{i}[/math]

с нечётными [math]i[/math] с помощью деления уравнения на [math]b^{(k)}_{i}[/math] выполняется его замена на уравнение

[math] (c^{(k)}_{i}/b^{(k)}_{i}) x^{(k)}_{i-1} + x^{(k)}_{i} + (a^{(k)}_{i}/b^{(k)}_{i}) x^{(k)}_{i+1} = f^{(k)}_{i}/b^{(k)}_{i}[/math]

Рисунок 1. Микрограф "узла" подготовительного шага прямого хода алгоритма циклической редукции

Уравнение

[math]b^{(k)}_{1} x^{(k)}_{1} + a^{(k)}_{1} x^{(k)}_{2} = f^{(k)}_{1}[/math]

аналогично меняется на уравнение

[math]x^{(k)}_{1} + (a^{(k)}_{1}/b^{(k)}_{1}) x^{(k)}_{2} = f^{(k)}_{1}/b^{(k)}_{1}[/math]:

а уравнение

[math]c^{(k)}_{n} x^{(k)}_{n-1} + b^{(k)}_{n} x^{(k)}_{n} = f^{(k)}_{n}[/math]

меняется на уравнение

[math](c^{(k)}_{n}/b^{(k)}_{n}) x^{(k)}_{n-1} + x^{(k)}_{n} = f^{(k)}_{n}/b^{(k)}_{n}[/math]:

Рисунок 2a. Микрограф "узла" с нечётным номером прямого хода алгоритма циклической редукции

Для каждого же из уравнений

[math]c^{(k)}_{i} x^{(k)}_{i-1} + b^{(k)}_{i} x^{(k)}_{i} + a^{(k)}_{i} x^{(k)}_{i+1} = f^{(k)}_{i}[/math]

с чётными [math]i[/math] (кроме [math]2[/math] и [math]n-2[/math]) выполняется, с учётом [math]x^{(k+1)}_{i/2} = x^{(k)}_{i}[/math] его замена на уравнение

[math]c^{(k+1)}_{i/2} x^{(k+1)}_{(i-2)/2} + b^{(k+1)}_{i/2} x^{(k+1)}_{i/2} + a^{(k+1)}_{i/2} x^{(k+1)}_{(i+2)/2} = f^{(k+1)}_{i/2}[/math]

Рисунок 2b. Микрограф "узла" с чётным номером прямого хода алгоритма циклической редукции


при этом

[math]c^{(k+1)}_{i/2} = - c^{(k)}_{i}(c^{(k)}_{i-1}/b^{(k)}_{i-1})[/math],

[math]a^{(k+1)}_{i/2} = - a^{(k)}_{i}(a^{(k)}_{i+1}/b^{(k)}_{i+1})[/math],

[math]b^{(k+1)}_{i/2} = b^{(k)}_{i} - c^{(k)}_{i}(a^{(k)}_{i-1}/b^{(k)}_{i-1}) - a^{(k)}_{i}(c^{(k)}_{i+1}/b^{(k)}_{i+1})[/math],

[math]f^{(k+1)}_{i/2} = f^{(k)}_{i} - c^{(k)}_{i}f^{(k)}_{i-1}/b^{(k)}_{i-1} - a^{(k)}_{i}f^{(k)}_{i-1}/b^{(k)}_{i-1}[/math].

Для 2го уравнения выполняется его замена на уравнение

[math]b^{(k+1)}_{1} x^{(k+1)}_{1} + a^{(k+1)}_{1} x^{(k+1)}_{(2} = f^{(k+1)}_{1}[/math]

при этом

[math]a^{(k+1)}_{1} = - a^{(k)}_{2}(a^{(k)}_{3}/b^{(k)}_{3})[/math],

[math]b^{(k+1)}_{1} = b^{(k)}_{2} - c^{(k)}_{2}(a^{(k)}_{1}/b^{(k)}_{1}) - a^{(k)}_{2}(c^{(k)}_{3}/b^{(k)}_{3})[/math],

[math]f^{(k+1)}_{1} = f^{(k)}_{2} - c^{(k)}_{2}f^{(k)}_{1}/b^{(k)}_{1} - a^{(k)}_{2}f^{(k)}_{1}/b^{(k)}_{1}[/math]


[math]n-1[/math]-е уравнение заменяется на

[math]c^{(k+1)}_{(n-1)/2} x^{(k+1)}_{(n-3)/2} + b^{(k+1)}_{(n-1)/2} x^{(k+1)}_{(n-1)/2} = f^{(k+1)}_{(n-1)/2}[/math]

при этом

[math]c^{(k+1)}_{(n-1)/2} = - c^{(k)}_{n-1}(c^{(k)}_{n-2}/b^{(k)}_{n-2})[/math],

[math]b^{(k+1)}_{(n-1)/2} = b^{(k)}_{n-1} - c^{(k)}_{n-1}(a^{(k)}_{n-2}/b^{(k)}_{n-2}) - a^{(k)}_{n-1}(c^{(k)}_{n}/b^{(k)}_{n})[/math],

[math]f^{(k+1)}_{(n-1)/2} = f^{(k)}_{n-1} - c^{(k)}_{n-1}f^{(k)}_{n-2}/b^{(k)}_{n-2} - a^{(k)}_{n-1}f^{(k)}_{n-2}/b^{(k)}_{n-2}[/math].

По окончании всех этих манипуляций размерность k+1-й СЛАУ оказывается равной [math](n-1)/2[/math].

Шаги повторяются до тех пор, пока после [math]m-1[/math] шагов редукции размерность СЛАУ не становится равной 1 и остаётся одно уравнение

[math]b^{(m-1)}_{1} x^{(m-1)}_{1} = f^{(m-1)}_{1}[/math]

1.2.2 Обратный ход редукции

Рисунок 3. Микрограф "узла" обратного хода алгоритма циклической редукции

Из последнего уравнения, полученного прямым ходом, вычисляется

[math]x^{(m-1)}_{1} = f^{(m-1)}_{1}/b^{(m-1)}_{1}[/math]

Теперь, последовательно уменьшая верхние индексы неизвестных, используется нечётные уравнения каждого шага для вычисления неизвестных с соотвествующими нечётными номерами. Чётные неизвестные получаются из тождеств [math]x^{(k)}_{i} = x^{(k+1)}_{i/2}[/math], а для нечётных [math]i[/math]

[math]x^{(k)}_{i} = f^{(k)}_{i}/b^{(k)}_{i} - (c^{(k)}_{i}/b^{(k)}_{i}) x^{(k)}_{i-1} - (a^{(k)}_{i}/b^{(k)}_{i}) x^{(k)}_{i+1} [/math]

c "левого края" системы будет

[math]x^{(k)}_{1} = f^{(k)}_{1}/b^{(k)}_{1} - (a^{(k)}_{1}/b^{(k)}_{1}) x^{(k)}_{2}[/math]

а с "правого"

[math]x^{(k)}_{n} = f^{(k)}_{n}/b^{(k)}_{n} - (c^{(k)}_{n}/b^{(k)}_{n}) x^{(k)}_{n-1}[/math]

После вычисления всех [math]x^{(0)}_{i}[/math] значения искомых неизвестных [math]x_{i} = x^{(k)}_{i}[/math].

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

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

[math] c^{(k)}_{i}/b^{(k)}_{i} , a^{(k)}_{i}/b^{(k)}_{i} , f^{(k)}_{i}/b^{(k)}_{i}[/math]

почти все используются для преобразований двух чётных уравнений, то при выделении "микровычислений", из которых следует составить шаги редукции и которые составляют его ядро, лучше отнести вычисления этих отношений к предыдущему шагу редукции. Таким образом, на "подготовительном шаге" микроядро будет для каждого нечётного [math]i[/math] состоять только из трёх делений (кроме [math]2[/math] и [math]n[/math] - там будет по 2 деления), а затем на каждом последующем шаге редукции для каждого нечётного [math]i[/math] - из двух умножений и двух вычислений выражений типа [math]a-bc-de[/math], с последующими тремя делениями (на "краях" часть этих операций отсутствует или урезана, но общую картину это не очень меняет). Для чётных [math]i[/math] деления в микроядре будут отсутствовать.

Что касается шагов обратного хода, то там для каждого [math]i[/math] рано или поздно выполняется одна операция типа [math]a-bc-de[/math] (на "краях" - типа [math]a-bc[/math]).

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

Рисунок 4. Граф алгоритма прямого хода циклической редукции при n=15. Светло-зелёным обозначены микроядра "подготовительного шага" с тремя (или меньше) делениями, синим - микроядра нечётных узов, сине-зелёным - ядра чётных узлов (без делений).
Рисунок 5. Граф алгоритма обратного хода циклической редукции при n=15. В вершинах - операции вида a-bc-de, в вершинах с 1 входящей дугой - вида a-bc. Показаны только дуги, передающие значения найденных неизвестных.

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

Обратный ход циклической редукции - уже практически полное "обратное" сдваивание, и на этом этапе разделение на независимые ветви может быть проделано после старта, если размножить необходимые данные.

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

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

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

Поскольку алгоритм циклической редукции обычно не предназначен для последовательного исполнения, то его "последовательную сложность" скорее следует считать только оценкой общего количества операций. Если взять только главный член, линейный по размеру, и опустить меньшие по порядку (логарифмические и константы), то получается, что в циклической редукции 3n делений, 8n умножений и 6n вычитаний/сложений. Таким образом, при классификации по последовательной сложности, алгоритм циклической редукции относится к алгоритмам с линейной сложностью.

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

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

У прямого хода циклической редукции макрограф представлен на Рис. 4, при этом графы разные макровершин представлены на Рис. 1, 2a, 2b. Как видно, после первого шага схему можно условно назвать "схемой страивания", поскольку связанные с каждым уравнением, остающимся в редуцированной СЛАУ, коэффициенты вычисляются по формулам из коэффициентов трёх соседних уравнений в нередуцированной СЛАУ.

Обратный ход имеет макрограф, представленный на Рис. 5, при этом графы макровершин представлены на Рис. 3 (крайние макровершины, имеющие только одну входящую дугу на макрографе, замещают недостающие данные нулями, соответственно уменьшая вычисления). Этот макрограф, в принципе, имеет макроструктуру, обратную к схеме сдваивания, поэтому её можно условно назвать "схемой раздваивания", с одной поправкой: все макровершины, кроме крайних, имеют по две входящих дуги, которые немного изменяют обычное "раздваивание".

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

Для выполнения циклической редукции в трёхдиагональной СЛАУ из [math]n[/math] уравнений с [math]n[/math] неизвестными в параллельном варианте требуется последовательно выполнить следующие ярусы:

  • [math]log_2 n[/math] ярусов делений, из них самый широкий ярус - первый, в нём [math]3n/2[/math] делений.
  • [math]2 log_2 n[/math] ярусов умножений и [math]4 log_2 n[/math]сложений/вычитаний.

Ширина ярусов экспоненциально, с показателем 2, убывает с их номерами, а потом, на обратном ходе, экспоненциально растёт, также с показателем 2. Эта неоднородность приводит к тому, что при практическом полном распараллеливании большую часть времени большинство устройств будет простаивать.

Таким образом, при классификации по высоте ЯПФ, прогонка относится к алгоритмам с логарифмической сложностью. При классификации по ширине ЯПФ её сложность будет линейной.

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

Входные данные: трёхдиагональная матрица [math]A[/math] (элементы [math]a_{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 Свойства алгоритма

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

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

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

Обычно циклическая редукция используется для решения СЛАУ с диагональным преобладанием. В этом случае гарантируется устойчивость алгоритма.

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

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

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

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

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

      subroutine cprogm (a,x,N,m) ! N=2^m-1
      real a(3,N), x(N)

      k=1
      k2=2
      kk=1  

      a(2,1)=1./a(2,1)
      a(2,N)=1./a(2,N)
      a(3,1)=a(3,1)*a(2,1)
      a(1,N)=a(1,N)*a(2,N)
      x(1) = x(1)*a(2,1)
      x(N) = x(N)*a(2,N)
      do i=3,N-2,2
         a(2,i) = 1./a(2,i)
         a(1,i) = a(1,i)*a(2,i)
         a(3,i) = a(3,i)*a(2,i)
         x(i) = x(i)*a(2,i)
      end do
      do j=1,m-2

         k=k2 ! k=2^j

         a(2,k) = a(2,k) - a(1,k)*a(3,k-kk)
         a(3,k) = -a(3,k)*a(3,k+kk)
         x(k) = x(k)- x(k-kk)*a(1,k-kk)
         a(2,N+1-k) = a(2,N+1-k) - a(1,N+1-k)*a(3,N+1-k-kk)
         a(1,N+1-k)= - a(1,N+1-k)*a(1,N+1-k-kk)
         x(N+1-k) = x(N+1-k)- x(N+1-k-kk)*a(1,N+1-k-kk)
         a(2,k) = a(2,k) - a(3,k)*a(1,k+kk)
         x(k) = x(k) - x(k+kk)*a(3,k+kk)
         a(2,N+1-k) = a(2,N+1-k) - a(3,N+1-k)*a(1,N+1-k+kk)
         x(N+1-k) = x(N+1-k)- x(N+1-k+kk)*a(3,N+1-k+kk)         

         k2=k2*2 ! k2=2^(j+1)

         do i = k2, N+1-k2, k
             a(2,i) = a(2,i) - a(1,i)*a(3,i-kk)
             x(i) = x(i)- x(i-kk)*a(1,i-kk)
             a(3,i) = -a(3,i)*a(3,i+kk)
             a(1,i)= - a(1,i)*a(1,i-kk)             
             a(2,i) = a(2,i) - a(3,i)*a(1,i+kk)
             x(i) = x(i)- x(i+kk)*a(3,i+kk)
         end do 

         a(2,k)=1./a(2,k)
         a(2,N+1-k)=1./a(2,N+1-k)
         a(3,k)=a(3,k)*a(2,k)
         a(1,N+1-k)=a(1,N+1-k)*a(2,N+1-k)
         x(k) = x(k)*a(2,k)
         x(N+1-k) = x(N+1-k)*a(2,N+1-k) 

         do i = k2+k, N+1-k2-k, k2
            a(2,i) = 1./a(2,i)
            a(1,i) = a(1,i)*a(2,i)
            a(3,i) = a(3,i)*a(2,i)
            x(i) = x(i)*a(2,i)
         end do
         
         kk=k  ! budet kk=2^(j-1)

      end do 

c        m-1 shag & start of reverse

      k=k2 ! k=2^(m-1) i kk=2^(m-2)

      a(2,k) = a(2,k) - a(1,k)*a(3,k-kk)
      x(k) = x(k)- x(k-kk)*a(1,k-kk)
      a(2,k) = a(2,k) - a(3,k)*a(1,k+kk)
      x(k) = x(k)- x(k+kk)*a(3,k+kk)
      a(2,k) = 1./a(2,k)
      x(k) = x(k)*a(2,k)

c           start reverse       

      x(kk) = x(kk) - a(3,kk)*x(k)

      x(k+kk) = x(k+kk) - a(1,k+kk)*x(k)

      do j=m-2, 1, -1

          k2 = k ! k2 = 2^(j+1)
          k = kk ! k = 2^j
          kk = kk/2 ! kk = 2^(j-1)

          x (kk) = x (kk) - a(3,kk)*x(kk+kk)
          x (N+1-kk) = x (N+1-kk) - a(1,N+1-kk)*x(N+1-k) 
          do i = k2+k, N+1-k-kk, k2
             x(i) = x(i) - a(1,i)*x(i-kk)  
          end do   
          do i = k2+k, N+1-k-kk, k2
             x(i) = x(i) - a(3,i)*x(i+kk) 
          end do  
      end do

      return
      end

Здесь необходимо отметить, что, несмотря на последовательность этой реализации, она вполне может быть исполнена и параллельно, нужно только как-то скомандовать имеющемуся в наличии компилятору, что циклы по i параллельны (например, спецкомментариями OpenMP).

2.2 Локальность данных и вычислений

По графу видно, что, хотя и налицо "местная" локальность по вычислениям, но от яруса к ярусу расстояние между запрашиваемыми данными растёт от первых и последних ярусов алгоритма к его середине, являющимся не только узким местом алгоритма, но и имеющим наибольшую длину передачи данных между устройствами. Поэтому, хотя в циклической редукции нет той избыточной локальности, которая притормаживает работу прогонки, это преимущество, как показали замеры[5] (см. Рисунок 6), не способно компенсировать недостатки общей локальности графа циклической редукции.

Рисунок 6. Отношение времени исполнения циклической редукции и монотонных прогонок на ПК. По оси асбсцисс отложено [math]m[/math] (порядок СЛАУ равен [math]2^m - 1[/math]), по оси ординат - отношение времени, затраченного на исполнение циклической редукции, к времени, затраченному на выполнение монотонной повторной прогонки. Система Core i3 + Windows 8, сборка кода компилятором GNU Fortran. На других системах получена качественно примерно такая же картина, с различием по выходу за границы кэша. Между синей и красной линиями - примерная область взаимной компенсации неарифметических факторов. При малых размерах сказывается неиспользование прогонками суперскалярности процессора.


2.2.1 Локальность реализации алгоритма

2.2.1.1 Структура обращений в память и качественная оценка локальности
Рисунок 7. Полный метод циклической редукции. Общий профиль обращений в память

На рис. 7 представлен профиль обращений в память для реализации полного метода циклической редукции. Данный профиль достаточно интересно. Как нам известно из описания самого алгоритма, он состоит из двух частей – прямого и обратного хода. На общем профиле синяя вертикальная черта отделяет эти два этапа. Рассмотрим их в отдельности.

На первом этапе явно видны отдельные итерации (первые 3 итерации разделены желтыми линиями), каждая из которых состоит из 4 параллельно выполняемых переборов данных. Отметим, что на каждой итерации в рамках каждого перебора задействованы данные из одного и того же диапазона. Число обращений в память в каждой итерации различается, и это позволяет предположить, что характер переборов также может отличаться.

Рассмотрим фрагменты двух итераций более детально. На рис. 8 и 9 соответственно представлены фрагменты 1 и 2, обозначенные на рис. 7 зеленым. Данные фрагменты содержат по 1000 обращений. Можно увидеть, что в целом все переборы в двух фрагментах выполняются с регулярным и достаточно небольшим шагом по памяти, однако интенсивность обращений в каждом случае разная. Причем при смене итерации меняется и интенсивность обращений к одному и тому диапазону данных.

Рисунок 8. Профиль обращений, фрагмент 1
Рисунок 9. Профиль обращений, фрагмент 2

Рассмотрим еще более детально небольшие части данных фрагментов, чтобы дальше проанализировать структуру представленных переборов. На рис. 10 представлены небольшие локальные области, выделенные зеленым на рис. 8 и 9. Можно увидеть, что локальная структура обращений в целом отличается. В частности, в области, выделенной на рис. 9, шаг по памяти больше, а также выполняется меньше повторных обращений, что говорит о более низкой пространственной локальности.

Рисунок 10. Локальная структура фрагментов 2 и 3

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

Отдельно рассмотрим самый конец первого этапа (фрагмент 3 на рис. 7), который показан на рис. 11. Здесь локальность обращений в память самая низкая, поскольку шаг по памяти в переборах становится наибольшим. В центре рисунка можно заметить область особенно низкой локальности (наиболее разрозненные обращения). Хотя данный фрагмент небольшого размера, он может серьезно снижать общую производительность взаимодействия с памятью.

Рисунок 11. Профиль обращений, фрагмент 3

Второй этап также имеет регулярную структуру. Рассмотрим на рис. 12 небольшую область более детально (фрагмент 4 на рис. 7). Из данного рисунка видно, что здесь также выполняются последовательные переборы, локальная структура которых может отличаться. Как и в рассмотренном на первом этапе случае, это приводит к различиям по пространственной и временной локальности. Стоит отметить, что угол наклона переборов на втором этапе в целом больше, чем на первом этапе (см. рис 7), что говорит о большем шаге по памяти, и, как следствие, более низкой пространственной локальности, однако это различие может нивелироваться различиями в локальной структуре.

Рисунок 12. Профиль обращений, фрагмент 4

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

2.2.1.2 Количественная оценка локальности

Оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.

На рисунке 13 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что производительность работы с памятью в данном случае низка и лишь чуть выше значения daps для такой неэффективной с точки зрения работы с памятью программы как тест RandomAccess (rand). В целом это соотносится с выводами, которые были сделаны нами ранее


Рисунок 13. Сравнение значений оценки daps

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

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

2.4 Масштабируемость алгоритма и его реализации

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

2.4.1 Масштабируемость алгоритма

2.4.2 Масштабируемость реализации алгоритма

Проведём исследование масштабируемости параллельной реализации циклической редукции согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов-2" Суперкомпьютерного комплекса Московского университета.

Набор и границы значений изменяемых параметров запуска реализации алгоритма:

  • число процессоров [2 : 256] с шагом степени двойки;
  • размер матрицы [64 : 33554432] с шагом степени двойки.

В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:

  • минимальная эффективность реализации 3.89e-09%;
  • максимальная эффективность реализации 0.00163%.

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

Рисунок 8. Параллельная реализация циклической редукции. Изменение производительности в зависимости от числа процессоров и размера матрицы.
Рисунок 9. Параллельная реализация циклической редукции. Изменение эффективности в зависимости от числа процессоров и размера матрицы.

Исследованная параллельная реализация на языке C

2.5 Динамические характеристики и эффективность реализации алгоритма

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

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

2.7 Существующие реализации алгоритма

Обычно циклическую редукцию как в блочном, так и в точечном варианте реализуют матфизики, решающие задачи с задачами, содержащими, например, дискретизацию оператора Лапласа. Её, несмотря на популярность, редко включают в пакеты программ. В основном это связано с тем, что простота схемы циклической редукции остаётся только для определённых размерностей задач, а универсальные решатели на её основе никто не делает. Для других размерностей могут быть придуманы разные, в том числе и более быстрые варианты циклической редукции[6], но у них очень сложная схема, и каждый раз нужно выполнять анализ на устойчивость.

3 Литература

  1. ↑ Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977.
  2. ↑ Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, 1984.
  3. ↑ 3,0 3,1 3,2 Ильин В.П., Кузнецов Ю.И. Трехдиагональные матрицы и их приложения. М.: Наука. Глав-ная редакция физико-математической литературы, 1985г., 208 с. Ошибка цитирования Неверный тег <ref>: название «IK3d» определено несколько раз для различного содержимого
  4. ↑ 4,0 4,1 Фролов А.В., Антонов А.С., Воеводин Вл.В., Теплов А.М. Сопоставление разных методов решения одной задачи по методике проекта Algowiki // Параллельные вычислительные технологии (ПаВТ’2016): труды международной научной конференции (г. Архангельск, 28 марта – 1 апреля 2016 г.). Челябинск: Издательский центр ЮУрГУ, 2016. С. 347-360.
  5. ↑ Фролов Н.А., Фролов А.В. Экспериментальные исследования влияния степени локальности алгоритмов на их быстродействие на примере решения трёхдиагональных СЛАУ // Труды 59й научной конференции МФТИ (21–26 ноября 2016 г., гг. Москва-Долгопрудный).
  6. ↑ А.В.Фролов "О коэффициенте при логарифме в критическом пути графа циклической редукции" // Суперкомпьютерные дни в России: Труды международной конференции (26-27 сентября 2016 г., г. Москва). – М.: Изд-во МГУ, 2016. с. 307-313.



Источник — https://algowiki-project.org/w/ru/index.php?title=Полный_метод_циклической_редукции&oldid=21045