Метод "Разделяй и властвуй" Завольсков/Землянский

Авторы статьи: Завольсков Владислав, Землянский Роман, 614 группа

1 ЧАСТЬ. Свойства и структура алгоритмов

1.1 Общее описание алгоритма

Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы - это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матриц, начиная с порядка n, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 25. Для матриц меньшего порядка (или если нужны только собственные значения) происходит автоматический переход к QR-итерации.

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

Пусть

[math] L = \begin{bmatrix} a_{1} & b_{1}&&&&& \\ b_{1} & \ddots & \ddots \\ & \ddots & a_{m-1} & b_{m-1} \\ && b_{m-1} & a_{m} & b_{m} \\ &&& b_{m} & a_{m+1} & b_{m+1} \\ &&&& b_{m+1} & \ddots \\ &&&&&& \ddots & b_{n-1} \\ &&&&&& b_{n-1} & a_{n} \\ \end{bmatrix} = \begin{bmatrix} a_{1} & b_{1} &&&&& \\ b_{1} & \ddots & \ddots \\ & \ddots & a_{m-1} & b_{m-1} \\ && b_{m-1} & a_{m} - b_{m} \\ &&&& a_{m+1} - b_{m} & b_{m+1} \\ &&&& b_{m+1} & \ddots \\ &&&&&& \ddots & b_{n-1} \\ &&&&&& b_{n-1} & a_{n} \\ \end{bmatrix} + \begin{bmatrix} &&&&& \\ &&b_{m} & b_{m} \\ &&b_{m} & b_{m} \\ &&&&& \\ \end{bmatrix} = [/math]

[math] = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m} * \begin{bmatrix} 0 \\ \vdots \\ 0 \\ 1 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix} = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} [/math]

Предположим, что нам известны спектральные разложения матриц [math]T_{1}[/math] и [math] T_{2} [/math]: [math] T_{i} = Q_{i} \Lambda_{i} Q_{i}^{T} [/math]. В действительности, они будут рекурсивно вычисляться тем же самым алгоритмом. Установим связь между собственными значениями матрицы Т и собственными значениями матриц [math]T_{1}[/math] и [math]T_{2}[/math]. Имеем:

[math] T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} = \begin{bmatrix} Q_{1} \Lambda_{1} Q_{1}^{T} & 0 \\ 0 & Q_{2} L_{2} Q_{2}^{T}\end{bmatrix} + b_{m}vv^{T} = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}(\begin{bmatrix} \Lambda_{1} & \\ & \Lambda_{2}\end{bmatrix} + b_{m}vv^{T})\begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix} [/math],

где

[math] u = \begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}v [/math]

так как [math]v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T[/math], получим матрицу, состоящую из последнего столбца матрицы [math] Q_{1}^{T}[/math] и первого столбца матрицы [math] Q_{2}^{T}[/math].

Следовательно, [math]T[/math] имеет те же собственные значения, что и пдобная ей матрица [math]D + \rho uu^{T}[/math], где [math]D = \begin{bmatrix} L_{1} & 0 \\ 0 & L_{2}\end{bmatrix}[/math] - диагональная матрица, [math]\rho = b_{m}[/math] - число, а [math]u[/math] - вектор.

Будем предполагать, не ограничивая общности, что диагональные элементы [math]d_{1}, \ldots, d_{n}[/math] матрицы [math]D[/math] упорядочены по убыванию: [math]d_{n} \lt = \ldots \lt =d_{1}[/math].

Чтобы найти собственные значения матрицы [math]D + \rho uu^{T}[/math], вычислим её характеристический многочлен, считая пока матрицу [math]D - \lambda I[/math] невырожденной. Тогда

[math]det(D + \rho uu^{T} - \lambda I) = det((D - \lambda I)(I + \rho (D- \lambda I)^{-1} uu^{T}))[/math].

Поскольку [math]D - \lambda I[/math] невырожденна, [math]det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 0[/math] тогда и только тогда, когда [math]\lambda[/math] - собственное значение. Заметим, что матрица [math]I + \rho (D - \lambda I)^{-1}uu^{T}[/math] получается из единичной добавлением матрицы ранга 1. Определитель такой матрицы легко вычислить.

Лемма 1. Справедливо равенство [math]det(I + xy^{T}) = 1 + y^{T}x[/math], где [math]x[/math] и [math]y[/math] - векторы.

Таким образом,

[math]det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 1 + \rho u^{T}(D - \lambda I)^{-1}u[/math]

[math] = 1 + \rho \sum_{i=1, n} \frac{u_{i}^{2}} {d_{i}-\lambda} = f(\lambda)[/math] ,

т.е. собственные значения матрицы [math]T[/math] есть корни так называемого векового уравнения [math]f(\lambda) = 0[/math]. Если все числа [math]d_{i}[/math] различны и все [math]u_{i} \lt \gt 0[/math] (случай общего положения), то [math]f(\lambda)[/math] имеет график типа показанного на рис. 1 (где [math]n = 4[/math] и [math]\rho \gt 0[/math]).

Можно видеть, что прямая [math]y = 1[/math] является горизонтальной асимптотой для этого графика, а прямые [math]\lambda = d_{i}[/math] есть вертикальные асимптоты. Поскольку (*ФОРМУЛА*), функция возрастает всюду, кроме точек [math]\lambda = d_{i}[/math]. Поэтому корни функции разделяются числами [math]d_{i}[/math] и ещё один корень находится справа от точки [math]d_{1}[/math] (на рис. 1 [math]d_{1} = 4[/math]). (При [math]\rho\lt 0[/math] функция [math]f(\lambda)[/math] всюду убывает и соответствующий корень находится слева от точки [math]d_{n}[/math]). Для функции [math]f(\lambda)[/math], монотонной и гладкой на каждом из интервалов [math](d_{i+1},d_{i})[/math], можно найти вариант метода Ньютона, который быстро и монотонно сходится к каждому из корней, если начальная точка взята в [math](d_{i+1},d_{i})[/math]. Нам достаточно знать, что на практике метод сходится к каждому собственному значению за ограниченное число шагов. Поскольку вычисление [math]f(\lambda)[/math] и [math]f^{'}(\lambda)[/math] стоит [math]O(n)[/math] флопов, для вычисления одного собственного значения достаточно [math]O(n)[/math], а для вычисления всех [math]n[/math] собственных значений матрицы [math]D + \rho uu^{T}[/math] требуется [math]O(n^{2})[/math] флопов. Для собственных векторов матрицы [math]D + \rho uu^{T}[/math] мы легко можем получить явные выражения.

Лемма 2. Если [math]\alpha[/math] - собственное значение матрицы [math]D + \rho uu^{T}[/math], то соответствующий вектор равен [math](D - \alpha I)^{-1}u[/math]. Поскольку матрица [math]D - \alpha I[/math] диагональная, для вычисления такого вектора достаточно [math]O(n)[/math] флопов.

Доказательство.

[math](D + \rho uu^{T})[(D - \alpha I)^{-1}u] = (D - \alpha I + \alpha I + \rho uu^{T})(D - \alpha I)^{-1}u[/math]

[math]=u + \alpha (D - \alpha I)^{-1}u + u[\rho u^{T}(D - \alpha I)^{-1}u] [/math]

[math]=u + \alpha(D - \alpha I)^{-1}u - u[/math]

                                     поскольку [math] \rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0 [/math]

[math]=\alpha [(D - \alpha I)^{-1}u][/math], что и требовалось.

Для вычисления по этой простой формуле всех [math]n[/math] собственных векторов требуется [math]O(n^{2})[/math] флопов. К сожалению, формула не обеспечивает численной устойчивости, так как для двух очень близких значений [math]\alpha_{i}[/math] может давать неортогональные приближенные собственные векторы [math]u_{i}[/math]. Потребовалось целое десятилетие для того, чтобы найти устойчивую альтернативу исходному описанию алгоритма. Снова детали будут обсуждаться позднее в данном разделе.

Алгоритм является рекурсивным.

1.3 Wiki Part

As with most eigenvalue algorithms for Hermitian matrices, divide-and-conquer begins with a reduction to tridiagonal form. For an [math]m \times m[/math] matrix, the standard method for this, via Householder reflections, takes [math]\frac{4}{3}m^{3}[/math] flops, or [math]\frac{8}{3}m^{3}[/math] if eigenvectors are needed as well. There are other algorithms, such as the Arnoldi iteration, which may do better for certain classes of matrices; we will not consider this further here.

In certain cases, it is possible to deflate an eigenvalue problem into smaller problems. Consider a block diagonal matrix

[math]T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix}.[/math]

The eigenvalues and eigenvectors of [math]T[/math] are simply those of [math]T_{1}[/math] and [math]T_{2}[/math], and it will almost always be faster to solve these two smaller problems than to solve the original problem all at once. This technique can be used to improve the efficiency of many eigenvalue algorithms, but it has special significance to divide-and-conquer.

The divide part of the divide-and-conquer algorithm comes from the realization that a tridiagonal matrix is "almost" block diagonal.

The size of submatrix [math]T_{1}[/math] we will call [math]n \times n[/math], and then [math]T_{2}[/math] is [math](m - n) \times (m - n)[/math]. Note that the remark about [math]T[/math] being almost block diagonal is true regardless of how [math]n[/math] is chosen (i.e., there are many ways to so decompose the matrix). However, it makes sense, from an efficiency standpoint, to choose [math]n \approx m/2[/math].

We write [math]T[/math] as a block diagonal matrix, plus a rank-1 correction:

The only difference between [math]T_{1}[/math] and [math]\hat{T}_{1}[/math] is that the lower right entry [math]t_{nn}[/math] in [math]\hat{T}_{1}[/math] has been replaced with [math]t_{nn} - \beta[/math] and similarly, in [math]\hat{T}_{2}[/math] the top left entry [math]t_{n+1,n+1}[/math] has been replaced with [math]t_{n+1,n+1} - \beta[/math].

The remainder of the divide step is to solve for the eigenvalues (and if desired the eigenvectors) of [math]\hat{T}_{1}[/math] and [math]\hat{T}_{2}[/math], that is to find the diagonalizations [math]\hat{T}_{1} = Q_{1} D_{1} Q_{1}^{T}[/math] and [math]\hat{T}_{2} = Q_{2} D_{2} Q_{2}^{T}[/math]. This can be accomplished with recursive calls to the divide-and-conquer algorithm, although practical implementations often switch to the QR algorithm for small enough submatrices.

The conquer part of the algorithm is the unintuitive part. Given the diagonalizations of the submatrices, calculated above, how do we find the diagonalization of the original matrix?

First, define [math]z^{T} = (q_{1}^{T},q_{2}^{T})[/math], where [math]q_{1}^{T}[/math] is the last row of [math]Q_{1}[/math] and [math]q_{2}^{T}[/math] is the first row of [math]Q_{2}[/math]. It is now elementary to show that

[math]T = \begin{bmatrix} Q_{1} & \\ & Q_{2} \end{bmatrix} \left( \begin{bmatrix} D_{1} & \\ & D_{2} \end{bmatrix} + \beta z z^{T} \right) \begin{bmatrix} Q_{1}^{T} & \\ & Q_{2}^{T} \end{bmatrix}[/math]

The remaining task has been reduced to finding the eigenvalues of a diagonal matrix plus a rank-one correction. Before showing how to do this, let us simplify the notation. We are looking for the eigenvalues of the matrix [math]D + w w^{T}[/math], where [math]D[/math] is diagonal with distinct entries and [math]w[/math] is any vector with nonzero entries.

If w_i is zero, ([math]e_i[/math],d_i) is an eigenpair of [math]D + w w^{T}[/math] since [math](D + w w^{T})e_i = De_i = d_i e_i[/math].

If [math]\lambda[/math] is an eigenvalue, we have:

[math](D + w w^{T})q = \lambda q[/math]

where [math]q[/math] is the corresponding eigenvector. Now

[math](D - \lambda I)q + w(w^{T}q) = 0[/math]

[math]q + (D - \lambda I)^{-1} w(w^{T}q) = 0[/math]

[math]w^{T}q + w^{T}(D - \lambda I)^{-1} w(w^{T}q) = 0[/math]

Keep in mind that [math]w^{T}q[/math] is a nonzero scalar. Neither [math]w[/math] nor [math]q[/math] are zero. If [math]w^{T}q[/math] were to be zero, [math]q[/math] would be an eigenvector of [math]D[/math] by [math](D + w w^{T})q = \lambda q[/math]. If that were the case, [math]q[/math] would contain only one nonzero position since [math]D[/math] is distinct diagonal and thus the inner product [math]w^{T}q[/math] can not be zero after all. Therefore, we have:

[math]1 + w^{T}(D - \lambda I)^{-1} w = 0[/math]

or written as a scalar equation,

[math]1 + \sum_{j=1}^{m} \frac{w_{j}^{2}}{d_{j} - \lambda} = 0.[/math]

This equation is known as the secular equation. The problem has therefore been reduced to finding the roots of the rational function defined by the left-hand side of this equation.

All general eigenvalue algorithms must be iterative, and the divide-and-conquer algorithm is no different. Solving the nonlinear secular equation requires an iterative technique, such as the Newton–Raphson method. However, each root can be found in O(1) iterations, each of which requires [math]\Theta(m)[/math] flops (for an [math]m[/math]-degree rational function), making the cost of the iterative part of this algorithm [math]\Theta(m^{2})[/math].

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

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

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

Сюда добавим про дефляцию и решение векового уравнения

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

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

proc dc_eig[math](T,Q,\Lambda)...[/math] по входной матрице [math]T[/math] вычисляются выходные матрицы [math]Q[/math] и [math]\Lambda[/math], такие, что [math]T = Q\Lambda Q^{T}[/math]

Если [math]T[/math] - матрица размера [math]1 x 1[/math]

1. Присвоить выходным параметрам значения [math] Q = 1, \Lambda = T[/math]

Иначе

1. Представить [math]T[/math] в виде [math] T = \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T} [/math]

2. call dc_eig[math](T_{1},Q_{1},\Lambda_{1})[/math]

3. call dc_eig[math](T_{2},Q_{2},\Lambda_{2})[/math]

4. Построить [math]D+\rho uu^{T}[/math] по [math] Lambda_{1},Lambda_{2}, Q_{1}, Q_{2}[/math]

5. Найти матрицу собственных значений [math]\Lambda[/math]

6. Найти матрицу собственных векторов [math]Q^{'}[/math] для матрицы [math]D+\rho uu^{T}[/math]

7. Построить матрицу собственных векторов [math]Q[/math] для матрицы [math]T[/math] :

[math] Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}* Q^{'} [/math]

8. Присвоить выходным параметрам значения [math]Q[/math] и [math]\Lambda[/math]

endif

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

Проанализируем сложность алгоритма. Пусть [math]t(n)[/math] - число флопов при обработке матрицы размера [math] n x n[/math] процедурой dc_eig. Тогда

[math] t(n) = 2t(n/2) [/math] два рекурсивных обращения к dc_eig<(math>T_{i},Q_{i},\Lambda_{i})</math>

[math] +O(n^{2})[/math] вычисление собственных значений матрицы [math]D+\rho uu^{T}[/math]

[math] +O(n^{2})[/math] вычисление собственных векторов матрицы [math]D+\rho uu^{T}[/math]

[math] +c*n^{3}[/math] вычисление матрицы [math]Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}*Q^{'}[/math]

Если [math] Q_{1}, Q_{2}[/math] и [math]Q^{'}[/math] рассматриваются как плотные матрицы и используется стандартный алгоритм матричного умножения, то константа [math] c [/math] в последней строке равна 1. Таким образом, именно это умножение составляет наиболее трудоёмкую часть алгоритма в целом. Игнорируя члены порядка [math]n^{2}[/math], получаем [math]t(n) = 2t(n/2) + cn^{3}[/math]. Решая это разностное уравнение, находим [math] t = c\frac{4}{3}n^{3} [/math]

На практике константа [math]c[/math] обычно гораздо меньше 1, потому что матрица [math]Q^{'}[/math] весьма разрежена вследствие явления, называемого дефляцией.

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

TBD

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

Здесь приводится оценка параллельной сложности алгоритма: числа шагов, за которое можно выполнить данный алгоритм в предположении доступности неограниченного числа необходимых процессоров (функциональных устройств, вычислительных узлов, ядер и т.п.). Параллельная сложность алгоритма понимается как высота канонической ярусно-параллельной формы [1]. Необходимо указать, в терминах каких операций дается оценка. Необходимо описать сбалансированность параллельных шагов по числу и типу операций, что определяется шириной ярусов канонической ярусно-параллельной формы и составом операций на ярусах.

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

Для описания ресурса параллелизма алгоритма (ресурса параллелизма информационного графа) необходимо указать ключевые параллельные ветви в терминах конечного и массового параллелизма. Далеко не всегда ресурс параллелизма выражается просто, например, через координатный параллелизм или, что то же самое, через независимость итераций некоторых циклов (да-да-да, циклы - это понятие, возникающее лишь на этапе реализации, но здесь все так связано… В данном случае, координатный параллелизм означает, что информационно независимые вершины лежат на гиперплоскостях, перпендикулярных одной из координатных осей). С этой точки зрения, не менее важен и ресурс скошенного параллелизма. В отличие от координатного параллелизма, скошенный параллелизм намного сложнее использовать на практике, но знать о нем необходимо, поскольку иногда других вариантов и не остается: нужно оценить потенциал алгоритма, и лишь после этого, взвесив все альтернативы, принимать решение о конкретной параллельной реализации. Хорошей иллюстрацией может служить алгоритм, структура которого показана на рис.2: координатного параллелизма нет, но есть параллелизм скошенный, использование которого снижает сложность алгоритма с [math]n\times m[/math] в последовательном случае до [math](n+m-1)[/math] в параллельном варианте.

Рассмотрим алгоритмы, последовательная сложность которых уже оценивалась в п.1.6. Параллельная сложность алгоритма суммирования элементов вектора сдваиванием равна [math]\log_2n[/math], причем число операций на каждом ярусе убывает с [math]n/2[/math] до [math]1[/math]. Параллельная сложность быстрого преобразования Фурье (базовый алгоритм Кули-Тьюки) для векторов с длиной, равной степени двойки - [math]\log_2n[/math]. Параллельная сложность базового алгоритма разложения Холецкого (точечный вариант для плотной симметричной и положительно-определенной матрицы) это [math]n[/math] шагов для вычислений квадратного корня, [math](n-1)[/math] шагов для операций деления и [math](n-1)[/math] шагов для операций умножения и сложения.

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

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

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

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

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

Вычислительная мощность алгоритма равна отношению числа операций к суммарному объему входных и выходных данных. Она показывает, сколько операций приходится на единицу переданных данных. Несмотря на простоту данного понятия, это значение исключительно полезно на практике: чем выше вычислительная мощность, тем меньше накладных расходов вызывает перемещение данных для их обработки, например, на сопроцессоре, ускорителе или другом узле кластера. Например, вычислительная мощность скалярного произведения двух векторов равна всего лишь [math]1[/math], а вычислительная мощность алгоритма умножения двух квадратных матриц равна [math]2n/3[/math].

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

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

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

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

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

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

"Длинные" дуги в информационном графе [1] говорят о потенциальных сложностях с размещением данных в иерархии памяти компьютера на этапе выполнения программы. С одной стороны, длина дуги зависит от выбора конкретной системы координат, в которой расположены вершины графа, а потому в другой системе координат они попросту могут исчезнуть (но не появится ли одновременно других длинных дуг?). А с другой стороны, вне зависимости от системы координат их присутствие может быть сигналом о необходимости длительного хранения данных на определенном уровне иерархии, что накладывает дополнительные ограничения на эффективность реализации алгоритма. Одной из причин возникновения длинных дуг являются рассылки скалярных величин по всем итерациям какого-либо цикла: в таком виде длинные дуги не вызывают каких-либо серьезных проблем на практике.

Для проектирования специализированных процессоров или реализации алгоритма на ПЛИС представляют интерес компактные укладки информационного графа [1], которые также имеет смысл привести в данном разделе.

2 ЧАСТЬ. Программная реализация алгоритма

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

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

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

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

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

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

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

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

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

На рис.3 и рис.4 показаны профили обращения в память для программ, реализующих разложение Холецкого и быстрое преобразование Фурье, по которым хорошо видна разница свойств локальности у данных алгоритмов.

Рис.3 Реализация метода Холецкого. Общий профиль обращений в память

Рис.4 Нерекурсивная реализация БПФ для степеней двойки. Общий профиль обращений в память

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

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

• представленный иерархически ресурс параллелизма, опирающийся на структуру циклических конструкций и на граф вызовов программы;
• комбинацию (иерархию) массового параллелизма и параллелизма конечного;
• возможные способы распределения операций между процессами/нитями;
• возможные способы распределения данных;
• оценку количества операций, объёма и числа пересылок данных (как общего числа, так и в пересчёте на каждый параллельный процесс);

и другие.

В этом же разделе должны быть даны рекомендации или сделаны комментарии относительно реализации алгоритма с помощью различных технологий параллельного программирования: MPI, OpenMP, CUDA или использования директив векторизации.

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

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

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

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

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

Рис.5 Масштабируемость классического алгоритма умножения плотных матриц в зависимости от числа процессоров и размера задачи

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

Это объемный раздел AlgoWiki, поскольку оценка эффективности реализации алгоритма требует комплексного подхода [5], предполагающего аккуратный анализ всех этапов от архитектуры компьютера до самого алгоритма. Основная задача данного раздела заключается в том, чтобы оценить степень эффективности параллельных программ, реализующих данный алгоритм на различных платформах, в зависимости от числа процессоров и размера задачи. Эффективность в данном разделе понимается широко: это и эффективность распараллеливания программы, это и эффективность реализации программ по отношению к пиковым показателям работы вычислительных систем.

Помимо собственно показателей эффективности, нужно описать и все основные причины, из-за которых эффективность работы параллельной программы на конкретной вычислительной платформе не удается сделать выше. Это не самая простая задача, поскольку на данный момент нет общепринятой методики и соответствующего инструментария, с помощью которых подобный анализ можно было бы провести. Требуется оценить и описать эффективность работы с памятью (особенности профиля взаимодействия программы с памятью), эффективность использования заложенного в алгоритм ресурса параллелизма, эффективность использования коммуникационной сети (особенности коммуникационного профиля), эффективность операций ввода/вывода и т.п. Иногда достаточно интегральных характеристик по работе программы, в некоторых случаях полезно показать данные мониторинга нижнего уровня, например, по загрузке процессора, кэш-промахам, интенсивности использования сети Infiniband и т.п. Хорошее представление о работе параллельной MPI-программы дают данные трассировки, полученные, например, с помощью системы Scalasca.

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

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

На практике это сделать можно по-разному: либо все свести в один текущий раздел, либо же соответствующие факты сразу включать в предшествующие разделы, где они обсуждаются и необходимы по смыслу. В некоторых случаях, имеет смысл делать отдельные варианты всей части II AlgoWiki применительно к отдельным классам архитектур, оставляя общей машинно-независимую часть I. В любом случае, важно указать и позитивные, и негативные факты по отношению к конкретным классам. Можно говорить о возможных вариантах оптимизации или даже о "трюках" в написании программ, ориентированных на целевые классы архитектур.

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

Для многих пар алгоритм+компьютер уже созданы хорошие реализации, которыми можно и нужно пользоваться на практике. Данный раздел предназначен для того, чтобы дать ссылки на основные существующие последовательные и параллельные реализации алгоритма, доступные для использования уже сейчас. Указывается, является ли реализация коммерческой или свободной, под какой лицензией распространяется, приводится местоположение дистрибутива и имеющихся описаний. Если есть информация об особенностях, достоинствах и/или недостатках различных реализаций, то это также нужно здесь указать. Хорошими примерами реализации многих алгоритмов являются MKL, ScaLAPACK, PETSc, FFTW, ATLAS, Magma и другие подобные библиотеки.

3 Литература

[1] Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. - СПб.: БХВ-Петербург, 2002. - 608 с.

[2] Воеводин В.В., Воеводин Вад.В. Спасительная локальность суперкомпьютеров //Открытые системы. - 2013. - № 9. - С. 12-15.

[3] Воеводин Вад.В., Швец П. Метод покрытий для оценки локальности использования данных в программах // Вестник УГАТУ. — 2014. — Т. 18, № 1(62). — С. 224–229.

[4] Антонов А.С., Теплов А.М. О практической сложности понятия масштабируемости параллельных программ// Высокопроизводительные параллельные вычисления на кластерных системах (HPC 2014): Материалы XIV Международной конференции -Пермь: Издательство ПНИПУ, 2014. С. 20-27.

[5] Никитенко Д.А. Комплексный анализ производительности суперкомпьютерных систем, основанный на данных системного мониторинга // Вычислительные методы и программирование. 2014. 15. 85–97.