Текущая версия на 13:28, 9 января 2017

Эта работа прошла предварительную проверку
Дата последней правки страницы: 09.01.2017
Данная работа соответствует формальным критериям.
Проверено IgorS.

Авторы статьи: Сергеев Никита (группа 614) (1.1, 1.2, 1.5, 1.6, 1.7, 1.8, 1.10, 2.4), Стрельцова Виктория (группа 614) (1.2, 1.3, 1.4, 1.6, 1.7, 1.9, 1.10, 2.4)

Метод "Разделяй и властвуй" вычисления собственных значений и векторов трёхдиагональной матрицы
Последовательный алгоритм
Последовательная сложность	$c\frac{4}{3}n^{3}$
Объём входных данных	$2n + 1$
Объём выходных данных	$n(n + 1)$
Параллельный алгоритм
Параллельная сложность	$O(n^{2.3})$ крайне редко: $O(n^{2})$

1 Свойства и структура алгоритма

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

Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы - является наиболее быстрым из существующих методов вычисления всех собственных значений и собственных векторов трехдиагональной матрицы, начиная с порядка n, который примерно равен 26. (Точное значение этого порогового порядка в конкретном случае зависит от компьютера.)

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

Для матриц порядка ниже 25 (или в тех случаях, когда требуются только собственные значения) автоматически происходит переход к методу QR-итерации.

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

Пусть

$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} = \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} \equiv \begin{bmatrix} T_{1} & 0 \\ 0 & T_{2}\end{bmatrix} + b_{m}vv^{T}$

В предположении, что известны спектральные разложения матриц $T_{1}$ и $T_{2}$ : $T_{i}=Q_{i} \Lambda_{i} Q_{i}^{T}$ . На самом деле, они будут вычисляться рекурсивно тем же самым алгоритмом.

Можно установить связь между собственными значениями матрицы Т и собственными значениями матриц $T_{1}$ и $T_{2}$ . Получается:

$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}$ ,

где $u = \begin{bmatrix} Q_{1}^{T} & 0 \\ 0 & Q_{2}^{T}\end{bmatrix}v$ - матрица, в первой строке которой находятся элементы последнего столбца матрицы $Q_{1}^{T}$ , а во второй - элементы первого столбца матрицы $Q_{2}^{T}$ , так как $v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T$ . Следовательно, $T$ имеет такие же собственные значения, что и подобная ей матрица $D + \rho uu^{T}$ , где $D = \begin{bmatrix} L_{1} & 0 \\ 0 & L_{2}\end{bmatrix}$ - диагональная матрица, $\rho = b_{m}$ - число, а $u$ - вектор. В дальнейшем предполагаем, не ограничивая общности, что диагональные элементы $d_{1}, \ldots, d_{n}$ матрицы $D$ упорядочены по не возрастанию: $d_{n} \lt = \ldots \lt =d_{1}$ .

Для того, чтобы найти собственные значения матрицы $D + \rho uu^{T}$ , потребуется вычислить её характеристический многочлен, считая пока матрицу $D - \lambda I$ невырожденной. Тогда получаем $det(D + \rho uu^{T} - \lambda I) = det((D - \lambda I)(I + \rho (D- \lambda I)^{-1} uu^{T}))$ .

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

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

Следовательно, получаем $det(I + \rho (D - \lambda I)^{-1}uu^{T}) = 1 + \rho u^{T}(D - \lambda I)^{-1}u$ $= 1 + \rho \sum_{i=1, n} \frac{u_{i}^{2}} {d_{i}-\lambda} \equiv f(\lambda)$ ,

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

Рис. 1. График функции

$f(\lambda) = 1 + \frac{0.5}{1 - \lambda} + \frac{0.5}{2 - \lambda} + \frac{0.5}{3 - \lambda} + \frac{0.5}{4 - \lambda}$

По рисунку можно заметить, что прямая $y = 1$ является горизонтальной асимптотой для данного графика, а прямые $\lambda = d_{i}$ есть вертикальные асимптоты. Так как $f^{'}(\lambda) = \rho \sum_{i=1, n} \frac{u_{i}^{2}} {(d_{i}-\lambda)^{2}}\gt 0$ , функция возрастает всюду, кроме точек $\lambda = d_{i}$ . Из этого следует, что корни функции разделяются числами $d_{i}$ и ещё один корень находится справа от точки $d_{1}$ (на рис. 1 $d_{1} = 4$ ). (При $\rho\lt 0$ функция $f(\lambda)$ всюду убывает и соответствующий корень будет находиться слева от точки $d_{n}$ ). Для функции $f(\lambda)$ , которая является монотонной и гладкой на каждом из интервалов $(d_{i+1},d_{i})$ , можно найти вариант метода Ньютона, который будет быстро и монотонно сходиться к каждому из корней, при условии, что начальная точка взята в $(d_{i+1},d_{i})$ . Нам достаточно знать тот факт, что на практике метод сходится к каждому собственному значению за строго ограниченное число шагов. Поскольку вычисление $f(\lambda)$ и $f^{'}(\lambda)$ стоит $O(n)$ флопов, для вычисления одного собственного значения достаточно $O(n)$ флопов, следовательно для вычисления всех $n$ собственных значений матрицы $D + \rho uu^{T}$ потребуется $O(n^{2})$ флопов. Для собственных векторов матрицы $D + \rho uu^{T}$ мы легко можем получить явные выражения.

 Лемма 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] флопов.

Доказательство. $(D + \rho uu^{T})[(D - \alpha I)^{-1}u] = (D - \alpha I + \alpha I + \rho uu^{T})(D - \alpha I)^{-1}u = u + \alpha (D - \alpha I)^{-1}u + u[\rho u^{T}(D - \alpha I)^{-1}u]=$

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

поскольку

$\rho u^{T}(D - \alpha I)^{-1}u + 1 = f(\alpha) = 0$ .

Что и требовалось доказать.

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

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

До сих пор полагалось, что все $d_{i}$ различны и все $u_{i}$ отличны от нуля. Если это будет не так, то вековое уравнение $f(\lambda)=0$ имеет $k$ вертикальных асимптот, где $k\lt n$ , а следовательно $k$ корней. Однако оказывается, что остальные $n - k$ собственных значений могут быть определены без каких-либо усилий: в случае, если $d_{i}=d_{i+1}$ или $u_{i}=0$ , легко показать, что $d_{i}$ является собственным значением и для матрицы $D + \rho uu^{T}$ . В такой ситуации мы говорим о дефляции. На практике выбирается некоторое пороговое значение и дефляция для числа $d_{i}$ регистрируется, в том случает, если в смысле этого порога $d_{i}$ достаточно близко к $d_{i+1}$ либо $u_{i}$ достаточно мало.

Основной выигрыш от использования дефляции состоит не в том, что убыстряется решение векового уравнения - этот этап в любом случае стоит лишь $O(n^{2})$ операций. Выигрыш заключается в ускорении матричного умножения на последнем шаге алгоритма. Действительно, если $u_{i}=0$ , то соответствующий собственный вектор есть i-й столбец $e_{i}$ единичной матрицы. Это означает, что $e_{i}$ является i-м столбцом в матрице $Q_{'}$ , следовательно при формировании матрицы $Q$ применяя левое умножение $Q_{1}$ на $Q_{2}$ вычисление i-го столбца не будет требовать никаких затрат. Аналогичное упрощение возможно в случае $d_{i} = d_{i+1}$ . При дефляции многих собственных значений устраняется большая часть работы, которая связана с матричным умножением.

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

Предположим, что некоторое $u_{i}$ , хотя и мало, все же недостаточно мало для того, чтобы была зарегистрирована дефляция. В этом случае применение метода Ньютона к решению векового уравнения встречается с затруднениями. Вспомним, что пересчет приближённого решения $u_{j}$ уравнения $f(\lambda) = 0$ в методе Ньютона основан на следующих положениях:

1. Вблизи точки $\lambda = \lambda_{j}$ функция $f(\lambda)$ аппроксимируется линейной функцией $l(\lambda)$ ; график есть прямая линия, касающаяся графика функции $f(\lambda)$ при $\lambda = \lambda_{j}$ .

2. В качестве $\lambda_{j+1}$ берётся нуль этого линейного приближения, т.е. $l(\lambda_{j+1})=0$ .

Функция, показанная на рис.1, не доставляет видимых трудностей методу Ньютона, поскольку вблизи каждого своего нуля $f(\lambda)$ достаточно хорошо аппроксимируется линейными функциями. Однако рассмотрим график функции на рис. 2. Она получена из функции на рис. 1 заменой значения .5 для $u_{i}^{2}$ на .001. Это новое значение недостаточно мало для того, чтобы вызвать дефляцию. График функции в левой части рис.2 визуально не отличим от её вертикальных и горизонтальных асимптот, поэтому в правой части укрупненно воспроизведён фрагмент графика, прилегающий к вертикальной асимптоте $\lambda = 2$ . Видно, что график слишком быстро "выполняет поворот" и для большей части значений $\lambda$ почти горизонтален. Поэтому, применяя метод Ньютона почти к любому начальному приближению $\lambda_{0}$ , мы получаем линейное приближение $l(\lambda)$ с почти горизонтальным графиком и малым положительным угловым коэффициентом. В результате $\lambda_{1}$ является отрицательным числом, огромным по абсолютной величине, которое совершенно бесполезно в качестве приближения к истинному корню.

Рис. 2. График функции

$f(\lambda) = 1 + \frac{10^{-3}}{1 - \lambda} + \frac{10^{-3}}{2 - \lambda} + \frac{10^{-3}}{3 - \lambda} + \frac{10^{-3}}{4 - \lambda}$

Чтобы найти выход из этого положения, можно модифицировать метод Ньютона следующим образом: раз $f(\lambda)$ нельзя хорошо приблизить линейной функцией $l(\lambda)$ , попробуем взять в качестве приближения какую-нибудь другую простую функцию $h(\lambda)$ . Нет ничего особого именно в прямых линиях: для метода Ньютона вместо $l(\lambda)$ можно взять любое приближение $h(\lambda)$ , значения и нули которого легко вычисляются. Функция $f(\lambda)$ имеет полюсы в точках $d_{i}$ и $d_{i+1}$ , которые определяют её поведение в соответствующих окрестностях. Поэтому при поиске корня в интервале $(d_{i+1}, d_{i})$ естественно выбрать функцию $h(\lambda)$ , также имеющую эти полюсы, т.е. функцию вида $h(\lambda)= \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3}$

Константы $c_{1},c_{2}$ и $c_{3}$ обеспечивающие, что $h(\lambda)$ есть приближение к $f(\lambda)$ , можно выбрать несколькими способами. Отметим, что если $c_{1},c_{2}$ и $c_{3}$ уже известны, то уравнение $h(\lambda)=0$ легко решается относительно $\lambda$ , поскольку сводится к эквивалентному квадратному уравнению $c_{1}(d_{i+1}-\lambda)+c_{2}(d_{i}-\lambda)+c_{3}(d_{i}-\lambda)(d_{i+1}-\lambda)=0$

Пусть $\lambda_{j}$ - приближённое значение корня. определим $c_{1},c_{2}$ и $c_{3}$ так, чтобы $\frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3} = h(\lambda) \approx f(\lambda) = 1 + \rho \sum_{k=1, n} \frac{u_{k}^{2}} {d_{k}-\lambda}$ для $\lambda$ в окрестности $\lambda_{j}$ . Заметим, что $f(\lambda) = 1 + \rho \sum_{k=1, i} \frac{u_{k}^{2}} {d_{k}-\lambda} + \rho \sum_{k=i+1, n} \frac{u_{k}^{2}} {d_{k}-\lambda} \equiv 1 + \psi_{1}(\lambda) + \psi_{2}(\lambda)$ .

Если $\lambda \in (d_{i+1},d_{i})$ , то $\psi_{1}(\lambda)$ есть сумма положительных слагаемых, а $\psi_{2}(\lambda)$ - сумма отрицательных. Поэтому и $\psi_{1}(\lambda)$ , и $\psi_{2}(\lambda)$ могут быть вычислены с высокой точностью; однако при их сложении вполне вероятно взаимное уничтожение верных разрядов и потеря относительной точности в сумме. Возьмем числа $c_{1}$ и $\hat{c_{1}}$ , такие, что функция $h_{1}(\lambda) \equiv \hat{c_{1}} + \frac{c_{1}}{d_{i}-\lambda}$ удовлетворяет условиям $h_{1}(\lambda_{j}) = \psi_{1}(\lambda_{j})$ и $h_{1}^{'}(\lambda_{j})=\psi_{1}^{'}(\lambda_{j})$ (*)

Это означает, что гипербола, являющаяся графиком функции $h_{1}(\lambda)$ , касается графика функции $\psi_{i}(\lambda)$ при $\lambda = \lambda_{j}$ . Два условия в (*) - это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола. Легко проверить, что $c_{1}=\psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})^{2}$ и $\hat{c_{1}}=\psi_{1}(\lambda_{j}) - \psi_{1}^{'}(\lambda_{j})(d_{i} - \lambda_{j})$ .

Подобным же образом выбираем $c_{2}$ и $\hat{c_{2}}$ так, чтобы функция $h_{2}(\lambda) \equiv \hat{c_{2}} + \frac{c_{2}}{d_{i+1}-\lambda}$ удовлетворяла условиям $h_{2}(\lambda_{j}) = \psi_{2}(\lambda_{j})$ и $h_{2}^{'}(\lambda_{j})=\psi_{2}^{'}(\lambda_{j})$

Наконец, полагаем $h(\lambda) = 1 + h_{1}(\lambda) + h_{2}(\lambda) = (1 + \hat{c_{1}} + \hat{c_{2}} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} \equiv c_{3} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda}$ .

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

Вычислительным ядром последовательной схемы решения будет являться вычисление матрицы $Q$ собственных векторов при помощи умножения матрицы $Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}$ на матрицу $Q^{'}$ . Эта операция имеет сложность $cn^{3}$ о чём будет говориться ниже, в разделе #Последовательная сложность алгоритма .

Данной операции предшествует вычисление собственных значений и векторов матрицы $D + \rho uu^{T}$ .

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

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

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

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

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

В случае, если $T$ - матрица размера $1$ x $1$ :

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

Иначе:

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

2. call dc_eig $(T_{1},Q_{1},\Lambda_{1})$

3. call dc_eig $(T_{2},Q_{2},\Lambda_{2})$

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

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

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

7. Построить матрицу собственных векторов $Q$ для матрицы $T$ : $Q = \begin{bmatrix} Q_{1} & 0 \\ 0 & Q_{2}\end{bmatrix}* Q^{'}$

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

endif

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

Рассмотрим сложность алгоритма. Пусть $t(n)$ - число флопов, которое потребуется при обработке матрицы размера $n * n$ процедурой 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]

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

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

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

На рисунке 3 изображена структура алгоритма.

Рис. 3. Дерево алгоритма "Разделяй и властвуй"

Операции:

$D$ - операция Divide разделения матрицы $T$ на матрицы $T_{1}$ и $T_{2}$

Индексы - размерность

$S_{k}$ - операция Conquer сбора матрицы $D$

$L_{k}$ - операция поиска собственных значений

$S^{'}$ - операция поиска собственных векторов матрицы $Q^{'}$

$S$ - операция поиска собственных векторов матрицы $Q$

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

В структуре алгоритма имеем два параллельных блока - вызов рекурсивных функций dc_eig для вычисления собственных значений и векторов матриц $T_{1}$ и $T_{2}$ - это единственная часть алгоритма, в которой возможно использовать параллелизм.

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

- В вершине дерева все процессы взаимодействуют.

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

- В листьях деревьев каждый процесс занимается вычислением независимо.

В случае реализации без параллелизма, функции dc_eig отрабатывают последовательно - сначала со входными параметрами $T_{1}, Q_{1}, \Lambda_{1}$ , затем - $T_{2}, Q_{2}, \Lambda_{2}$ .

Таким образом, благодаря параллелизму высота алгоритма параллельного составляет $log_{2}(n)$ , тогда как при последовательной реализации высота алгоритма равна сумме геометрической прогрессии $1+2+4+8+...+2^{n-1} = 2^{n} - 1$ , что значительно сокращает число итераций алгоритма.

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

При размерности матрицы больше 3 параллельная реализация позволяет сократить количество этапов работы алгоритма больше, чем в два раза.

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

Входные данные: трёхдиагональная симметрическая матрица, описанная в разделе #Математическое описание алгоритма.

$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}$

Объём входных данных: $2n + 1$ (в силу симметричности достаточно хранить только диагональ и над/под диагональные элементы).

Выходные данные: Вектор собственных значений длины $n$ и $n$ собственных векторов длины $n$ , соответствующих собственным значениям исходной матрицы.

Объём выходных данных: $n(n + 1)$

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

Данный алгоритм является самым быстрым алгоритмом среди существующих: Разделяй и властвуй, QR, Бисекция и обратная итерация.

Последовательная сложность данного алгоритма: $t = c\frac{4}{3}n^{3}$

Она удовлетворяет соотношению $t(n) = n^3 +2t(\frac{n}{2})$

Параллельная сложность алгоритма (за счет дефляции): $O(n^{2.3})$ ; в крайне редких случаях: $O(n^{2})$ .

Особенности алгоритма:

1. Дефляция (#Дефляция)

2. Адаптированный метод Ньютона (#Макроструктура алгоритма)

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

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

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

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

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

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

Тестирование реализации алгоритма проводилось на суперкомпьютере "Ломоносов" с использованием 1, 2, 4, 8, 16 процессоров на матрицах размера от 4 до 4096. Реализация произведена на библиотеке OpenMP.

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

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

Рис. 4. Зависимость времени работы от размера матрицы

Параметры компиляции и запуска:

ssh compiler 
module add slurm
module add openmpi/1.5.5-icc
g++ _scratch/DivideAndConquer.c++ -fopenmp -std=c++0x -o DC 
sbatch -p test -n1 run ./DC

Код реализации: https://drive.google.com/open?id=0B0G2UV7cuLwHQk9xd1lyeDNJOXc

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

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

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

Единственный найденный пример реализации описан в Francoise Tisseury, Jack Dongarra, "Parallelizing the Divide and Conquer Algorithm for the Symmetric Tridiagonal Eigenvalue Problem on Distributed Memory Architectures" ([3] с. 7, 19-20).

3 Литература

[1] Дж. Деммель, «Вычислительная линейная алгебра» //С. 230-235

[2] Doron Grill and Eitan Tadmor AN $O(N2)$ METHOD FOR COMPUTING THE EIGENSYSTEM OF $N$ x $N$ SYMMETRIC TRIDIAGONAL MATRICES BY THE DIVIDE AND CONQUER APPROACH

[3] Francoise Tisseury, Jack Dongarra, A Parallel Divide and Conquer Algorithm for the Symmetric Eigenvalue Problem on Distributed Memory Architectures

[4] Francoise Tisseury, Jack Dongarra, Parallelizing the Divide and Conquer Algorithm for the Symmetric Tridiagonal Eigenvalue Problem on Distributed Memory Architectures

[5] Алгоритм "Разделяй и властвуй" - Wikipedia

Участник:Vika3410/Метод Разделяй и властвуй: различия между версиями