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

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

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

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

Метод разделяй и властвуй вычисления собственных значений и векторов трёхдиагональной матрицы - это наиболее быстрый из существующих методов, если нужны все собственные значения и собственные векторы трехдиагональной матриц, начиная с порядка n, примерно равного 26. (Точное значение этого порогового порядка зависит от компьютера.) Его численно устойчивая реализация весьма не тривиальна. В самом деле, хотя впервые метод был предложен еще в 1981 г., "правильный" способ его реализации был найден лишь в 1992 г. Этот способ воплощен LAPACK-программами ssyevd (для плотных матриц) и sstevd (для трехдиагональных матриц). В них стратегия "разделяй-и-влавствуй" используется для матриц порядка, большего чем 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}$. Имеем:

 [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]

так как $v = \begin{bmatrix} 0 , \ldots , 0 , 1 , 1 , 0 \ldots , 0 \end{bmatrix}^T$, получим матрицу, состоящую из последнего столбца матрицы $Q_{1}^{T}$ и первого столбца матрицы $Q_{2}^{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$ невырожденной. Тогда

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

Поскольку $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] - векторы.

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

 [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} \equiv f(\lambda)[/math] , 

т.е. собственные значения матрицы $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] флопов.

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

 [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], что и требовалось.

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

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

1.2.1 Дефляция

До сих пор полагалось, что все $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}$. При дефляции многих собственных значений устраняется большая часть работы, связанной с матричным умножением.

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

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

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

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

1.4.1 Решение векового уравнения

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

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

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

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

Функция, показанная на рис.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)$, также имеющую эти полюсы, т.е. функцию вида

[math]h(\lambda)= \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda} + c_{3}[/math]

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

[math]c_{1}(d_{i+1}-\lambda)+c_{2}(d_{i}-\lambda)+c_{3}(d_{i}-\lambda)(d_{i+1}-\lambda)=0[/math]

Пусть $\lambda_{j}$ - приближённое значение корня. определим $c_{1},c_{2}$ и $c_{3}$ так, чтобы

[math]\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} [/math]

для $\lambda$ в окрестности $\lambda_{j}$. Заметим, что

[math]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)[/math].

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

[math]h_{1}(\lambda) \equiv \hat{c_{1}} + \frac{c_{1}}{d_{i}-\lambda}[/math]   удовлетворяет условиям  [math]h_{1}(\lambda_{j}) = \psi_{1}(\lambda_{j})[/math] и [math]h_{1}^{'}(\lambda_{j})=\psi_{1}^{'}(\lambda_{j})[/math] (1)

Это означает, что гипербола, являющаяся графиком функции $h_{1}(\lambda)$, касается графика функции $\psi_{i}(\lambda)$ при $\lambda = \lambda_{j}$. Два условия в (1) - это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола. Легко проверить, что $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}}$ так, чтобы функция

[math]h_{2}(\lambda) \equiv \hat{c_{2}} + \frac{c_{2}}{d_{i+1}-\lambda}[/math]

Удовлетворяла условиям

 [math]h_{2}(\lambda_{j}) = \psi_{2}(\lambda_{j})[/math] и [math]h_{2}^{'}(\lambda_{j})=\psi_{2}^{'}(\lambda_{j})[/math] (2)

Наконец, полагаем

[math]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} [/math]

[math]\equiv c_{3} + \frac{c_{1}}{d_{i}-\lambda} + \frac{c_{2}}{d_{i+1}-\lambda}[/math].

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

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

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

Если $T$ - матрица размера $1$ x $1$

 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.6 Последовательная сложность алгоритма

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

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

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

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

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

Если $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 Информационный граф

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

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

Рис. 4. Детальное описание одного блока алгоритма

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

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

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

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

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

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

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

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

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

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

1. Использование дефляции (#Дефляция)

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

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

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

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

TBD

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

Большого количества реализаций данного алгоритма найдено не было: единственный пример описан на страницах 7, 19 и 20 в этой статье, которая указана в списке литературы.

3 Литература

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

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

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

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

[5] 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