Алгоритм Ланцоша для точной арифметики (без переортогонализации)
Основные авторы статьи (разделы 1, 2.4, 2.7):
Башлыков Л.А. и
Волков Н.И.
Авторы отдельных разделов (2.4, 2.7):
Григорьев М.А.,
Заспа А.Ю.,
Левин А.Д.,
Слюсарь Д.Р. и
Фролов А.А.
В работе также принимали участие:
Антонов А.С. и
Коньшин И.Н.
| Алгоритм Ланцоша для точной арифметики (без переортогонализации) | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]O(n^2) + O(k^3)[/math] - на одну итерацию и [math]O(kn^2)[/math] - суммарно для [math]k[/math] итераций |
| Объём входных данных | [math]n^2 + n + 1[/math], при экономии памяти затраты сокращаются до [math]\frac{n(n - 1)}{2} + n + 1[/math] |
| Объём выходных данных | [math]k[/math], а если нужны собственные векторы, то [math]k + kn[/math] |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]O(n)[/math] |
| Ширина ярусно-параллельной формы | [math]O(n^2)[/math] |
Содержание
- 1 Свойства и структура алгоритма
- 1.1 Общее описание алгоритма
- 1.2 Математическое описание алгоритма
- 1.3 Вычислительное ядро алгоритма
- 1.4 Макроструктура алгоритма
- 1.5 Схема последовательной реализации алгоритма
- 1.6 Последовательная сложность алгоритма
- 1.7 Информационный граф
- 1.8 Ресурс параллелизма алгоритма
- 1.9 Входные и выходные данные алгоритма
- 1.10 Свойства алгоритма
- 2 Программная реализация алгоритма
- 2.1 Особенности реализации последователього алгоритма
- 2.2 Локальность данных и вычислений
- 2.3 Возможные способы и особенности параллельной реализации алгоритма
- 2.4 Масштабируемость алгоритма и его реализации
- 2.5 Динамические характеристики и эффективность реализации алгоритма
- 2.6 Выводы для классов архитектур
- 2.7 Существующие реализации алгоритма
- 3 Литература
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Метод Ланцоша — процедура, применяемая для нахождения [math]k[/math] собственных значений симметричной матрицы произвольного размера путём вычисления их приближений. Суть метода заключается в итеративном построении специальной матрицы [math]T_{j}[/math], а затем вычислении её собственных значений и собственных векторов. Общее описание алгоритма можно найти, например, в книге [1].
Здесь мы рассмотрим сильно [math]x + y[/math] упрощенный вариант алгоритма Ланцоша, который подразумевает отсутствие влияния ошибок округления на вычислительный процесс. На практике используются модификации алгоритма Ланцоша, обладающие устойчивостью и без этого допущения ( например, алгоритм Ланцоша с полной переортогонализацией ).
1.2 Математическое описание алгоритма
Математическое описание алгоритма Ланцоща для вычисления собственных значений симметричной матрицы А требует 1) математического описания метода Ланцоща для построения крыловского подпространства и 2) математического описания т.н. процедуры Рэлея-Ритца.
1.2.1 Метод Ланцоща для построения крыловского подпространства
Пусть дано уравнение [math]Ax = b[/math], причём вектор [math]b[/math] известен и имеется
способ вычисления [math]Ax[/math].
Вычислим [math]y_{1} = b[/math] , [math]y_{2} = Ab[/math], [math]y_{3} = A^{2}b[/math],...,[math]y_{n} = A^{(n - 1)}b[/math], где [math]n[/math] -
порядок матрицы [math]A[/math]. Определим матрицу [math]K = [y_{1},...,y_{n}][/math].
Тогда [math]AK = [Ay_{1},...,Ay_{n - 1}, Ay_{n}] = [y_{2},...,y_{n},A^n{y_{1}}][/math].
В предположении, что матрица [math]K[/math] - невырождена, можно вычислить вектор [math]c = -K^{-1}A^{ny_{1}}[/math].
Поскольку первые [math]n - 1[/math] столбцов в [math]AK[/math] совпадают с последними [math]n - 1[/math] столбцами в
[math]K[/math], то справедливо: [math]AK = K[e_{2},e_{3},...,e_{n},-c] = KC[/math], то есть [math]K^{-1}AK = C[/math], где [math]C[/math] - верхняя хессенбергова матрица.
Заменим [math]K[/math] ортогональной матрицей [math]Q[/math] так, что при любом [math]k[/math] линейные оболочки первых [math]k[/math] столбцов в [math]K[/math] и [math]Q[/math] - одно и то же подпространство. Это подпространство называется крыловским и точно определяется как [math]\kappa_{k}(A, b)[/math] - линейная оболочка векторов [math]b[/math],[math]Ab[/math],...,[math]A^{(k - 1)}b[/math]. Суть алгоритма Ланцоша заключается в вычислении стольких первых столбцов в [math]Q[/math], сколько необходимо для получения требуемого приближения к решению [math]Ax = b [/math] ([math]Ax = {\lambda}x[/math]).
Используем QR-разложение матрицы [math]K = QQ[/math]
Тогда [math]K^{-1}AK = (R^{-1}Q^{T})A(QR) = C[/math], откуда [math]Q^{T}AQ = RCR^{-1} = H[/math]
Матрица [math]H[/math] является верхней хессенберговой в силу того, что верхней хессенберговой является матрица [math]C[/math], а матрицы [math]R[/math] и [math]R^{-1}[/math] - верхние треугольные.
Однако алгоритм Ланцоша подразумевает симметричность матрицы [math]A[/math]. Положим [math]Q = [Q_{k}, Q_{u}][/math], где [math]Q_{k} = [q_{1},...,q_{k}][/math] и [math]Q_{u} = [q_{k+1},...,q_{n}][/math]. В итоге получается, что [math]H = Q^{T}AQ = [Q_{k},Q_{u}]^{T}A[Q_{k},Q_{u}][/math] Из симметричности [math]A[/math], таким образом, следует симметричность и трехдиагональность матрицы [math]H = T[/math]. Теперь можно вывести две основные формулы, используемые в методе Ланцоша:
- [math]AQ_{j} = \beta_{j - 1}q_{j - 1} + \alpha_{j}q_{j} + \beta_{j}q_{j + 1}[/math] - получается приравниванием столбцов [math]j[/math] в обеих частях равенства [math]AQ = QT[/math].
- [math]q_{j}^{T}Aq_{j} = \alpha_{j}[/math] - получается домножением обеих частей предыдущего соотношения на [math]q_{j}^{T}[/math] и учетом ортогональности [math]Q[/math].
1.2.2 Процедура Рэлея-Ритца
Пусть [math]Q = [Q_{k}, Q_{u}][/math] - произвольная ортогональная матрица порядка [math]n[/math], причём [math]Q_{k}[/math] и [math]Q_{u}[/math] имеют размеры [math]nk[/math] и [math]n(n - k)[/math].
Пусть [math]T = Q^{T}AQ = [Q_{k},Q_{u}]^{T}A[Q_{k},Q_{u}][/math]. Обратим внимание, что в случае [math]k = 1[/math] выражение [math]T_{k} = Q_{k}^{T}AQ_{k}[/math] -
отношение Рэлея [math]\rho(Q_{1},A)[/math], то есть число [math](Q_{1}^{T}AQ_{1}) / (Q_{1}^{T}Q_{1})[/math]. Определим теперь суть
процедуры Рэлея-Ритца:
- Процедура Рэлея-Ритца - интерпретация собственных значений матрицы [math]T = Q^{T}AQ[/math] как приближений к собственным значениям матрицы [math]A[/math]. Эти приближения называются числами Ритца. Пусть [math]T_{k} = V{\lambda}V^{T}[/math] - спектральное разложение матрицы [math]T_{k}[/math]. Столбцы матрицы [math]Q_{k}V[/math] рассматриваются как приближения к соответствующим собственным векторам и называются векторами Ритца.
Существует несколько важных фактов о процедуре Рэлея-Ритца, характеризующих получаемые приближения собственных значений матрицы [math]A[/math].
- Минимум величины [math]||AQ_{k} - Q_{k}R||_{2}[/math] по всем симметричным [math]k^2[/math] матрицам [math]R[/math] достигается при [math]R = T_{k}[/math]. При этом [math]||AQ_{k} - Q_{k}R||_{2} = ||T_{ku}||_{2}[/math]. Пусть [math]T_{k} = V{\lambda}V^{T}[/math] - спектральное разложение матрицы [math]T_{k}[/math]. Минимум величины [math]||AP_{k} - P_{k}D||_{2}[/math], когда [math]P_{k}[/math] пробегает множество [math]nk [/math] матрицы с ортонормированными столбцами, таких, что [math]span(P_{k}) = span(Q_{k})[/math], а [math]D[/math] пробегает множество диагональных [math]k^2[/math] матриц также равен [math]||T_{ku}||_{2}[/math] и достигается для [math]P_{k} = Q_{k}V[/math] и [math]D = \lambda[/math].
Пусть [math]T_{k}[/math], [math]T_{ku}[/math], [math]Q_{k}[/math] - те же самые матрицы. Пусть [math]T_{k} = V{\lambda}V^{T}[/math] - всё то же спектральное разложение, причём [math]V = [v_{1},...,v_{k}][/math] - ортогональная матрица, а [math]\lambda = diag(\theta_{1},...,\theta_{k})[/math]. Тогда:
- Найдутся [math]k[/math] необязательно наибольших собственных значений [math]\alpha_{1},...,\alpha_{k}[/math] матрицы [math]A[/math], такие, что [math]|\theta_{i} - \alpha_{i}| \lt = ||T_{ku}||_{2}[/math] для [math]i = 1,...,k[/math]. Если [math]Q_{k}[/math] вычислена методом Ланцоша,то правая часть этого выражения равняется [math]\beta_{k}[/math], то есть единственному элементу блока [math]T_{ku}[/math], который может быть отличен от нуля. Указанный эдемент находится в правом верхнем углу блока.
- [math]||A(Q_{k}v_{i}) - (Q_{k}v_{i})\theta_{i}||_{2} = ||T_{ku}v_{i}||_{2}[/math]. Таким образом, разность между числом Ритца [math]\theta_{j}[/math] и некоторым собственным значением [math]\alpha[/math] матрицы [math]A[/math] не превышает величины [math]||T_{ku}v_{i}||_{2}[/math], которая может быть много меньше, чем [math]||T_{ku}||_{2}[/math]. Если матрица [math]Q_{k}[/math] вычислена алгоритмом Ланцоша, то эта величина равна [math]\beta_{k}|v_{i}(k)|[/math], где [math]v_{i}(k)[/math] - k-ая компонента вектора [math]v_{i}[/math]. Таким образом, можно дещево вычислить норму невязки - левой части выражения, не выполняя ни одного умножения вектора на матрицы [math]Q_{k}[/math] или [math]A[/math].
- Не имея информации о спектре матрицы [math]T_{u}[/math], нельзя дать какую-либо содержательную оценку для погрешности в векторе Ритца [math]Q_{k}v_{i}[/math]. Если известно, что [math]\theta_{j}[/math] отделено расстоянием, не меньшим [math]g[/math], от прочих собственных значений матриц [math]T_{k}[/math], [math]T_{u}[/math] и матрица [math]Q_{k}[/math] вычислена алгоритмом Ланцоша, то угол [math]\theta[/math] между [math]Q_{k}v_{i}[/math] и точным собственным вектором [math]A[/math] можно оценить формулой [math]\frac{1}{2}\sin{2\theta} \lt = \frac{b_{k}}{g}[/math]
В алгоритме Ланцоша из ортонормированных векторов строится матрица [math]Q_{k} = [q_{1},...,q_{k}][/math] и в качестве приближенных собственных значений [math]A[/math] принимаются числа Ритца - собственные значения симметричной трёхдиагональной матрицы [math]T_{k} = Q_{k}^{T}AQ_{k}[/math].
1.3 Вычислительное ядро алгоритма
Вычислительное ядро алгоритма Ланцоща состоит в получении на каждой итерации очередного промежуточного вектора [math]z[/math], получаемого путём умножения исходной матрицы [math]A[/math] на вектор Ланцоша [math]q_{j}[/math], полученный на предыдущей итерации.
Также при значениях [math]k[/math], сопоставимых с [math]n[/math], процедура вычисления собственных значений и собственных векторов симметричной трёхдиагональной матрицы по некоторому выбранному алгоритму подразумевает значительный объём расчётов и может считаться вычислительным ядром. Однако на практике метод используется прежде всего при малых [math]k[/math].
1.4 Макроструктура алгоритма
Метод Ланцоша соединяет метод Ланцоша для построения крыловского подпространства в процедурой Рэлея-Ритца, подразумевающей вычисление собственных значений симметричной трёхдиагональной матрицы. Первая часть алгоритма представляет собой строго последовательные итерации, на каждой из которых вначале строится очередной столбец матрицы [math]Q_{j}[/math] из первых [math]j[/math] ортонормированных векторов Ланцоша и матрица [math]T_{j} = Q^T_{j}AQ_{j}[/math] порядка [math]j[/math]. Итерационный процесс завершается, когда [math]j = k[/math]. Вторая часть алгоритма представляет собой поиск собственных значений и собственных векторов симметричной трёхдиагональной матрицы [math]T_{j}[/math], что реализуется отдельным алгоритмом, на практике используется QR-алгоритм.
Особенность методов крыловского подпространства состоит в предположении, что матрица [math]A[/math] доступна в виде "черного ящика", а именно подпрограммы, вычисляющей по заданному вектору [math]z[/math] произведение [math]y = Az[/math] (а также, возможно, произведение [math]y = A^{T}z[/math], если матрица [math]A[/math] несимметрична). Другими словами, прямой доступ к элементам матрицы и их изменение не используются, т.к.
- Умножение матрицы на вектор - наиболее дешевая нетривиальная операция, которую можно проделать с разреженной матрицей: в случае разреженной [math]A[/math], содержащей [math]m[/math] ненулевых элементов, для матрично-векторного умножения нужны [math]m[/math] скалярных умножений и не более [math]m[/math] сложений.
- [math]A[/math] может быть не представлена в виде матрицы явно, а доступна именно через подпрограмму для вычисления произведений [math]Ax[/math].
Таким образом, как QR-алгоритм поиска собственных значений, собственных векторов матрицы [math]T_{j}[/math] и оценки погрешности в них результирующей трехдиагональной матрицы, так и умножение матрицы на вектор внутри каждой итерации могут быть рассмотрены в качестве отдельных макроопераций. Проблемой является то, что именно эти манипуляции порождают основную вычислительную сложность алгоритма. Поскольку матрица, к которой применяется QR-алгоритм, имеет порядок [math]k[/math] и на практике бывает мала, сделаем допущение о плотности матрицы [math]A[/math]. Тогда умножение этой матрицы на вектор внутри каждой итерации нужно будет реализовывать в рамках самого алгоритма. Еще одним важным замечанием является момент вычисления собственных значений матрицы [math]T_{j}[/math]. Её собственные значения, полученные на итерации [math]p[/math], никак не используются на итерации [math]p + 1[/math], поэтому целесообразно вынести эту процедуру за основной цикл. Однако на практике в силу малых размеров матрицы [math]T_{j}[/math] вычисление соответствующих величин производится непосредственно в теле основного цикла, что позволяет сразу оценить достигнутую точность вычислений.
1.5 Схема последовательной реализации алгоритма
Алгоритм итерационный, на каждой итерации выполняется вычисление очередного столбца матрицы [math]Q[/math], а также вычисление элементов очередной пары ( диагонального и околодиагонального ) элементов матрицы [math]T[/math] с использованием значений столбца матрицы [math]Q[/math], полученного на предыдущей итерации и значения околодиагонального элемента матрицы [math]T[/math] с предыдущей итерации. Вектор [math]z[/math] носит вспомогательный характер. Коэффициенты [math]\alpha_{j}[/math] соответствуют диагональным элементам матрицы [math]T_{j}[/math], коэффициенты [math]\beta_{j}[/math] – наддиагональным и поддиагональным элементам той же матрицы. Далее приводится пример псевдокода. По окончании итераций вычислеяются собственные значения и собственные векторы матрицы [math]T_{j}[/math], что в псевдокоде не отражается ( т.к. это отдельный алгоритм ).
Вычислить столбец q1 матрицы Qk q1 = b / ||b||
q0 - нулевой вектор
beta_0 = 0
foreach (j = 1 .. k )
{
z = A * qj
alpha_j = qj * z
z = z – alpha_j * qj – beta_(j-1) * q(j – 1)
beta_j = ||z||
if ( beta_j == 0 )
{
break // Все ненулевые собственные значения найдены. Выход может быть досрочным.
}
else
{
q(j + 1) = z/beta_j
}
}
Вычислить собственные значения и собственные векторы симметричной трехдиагональной матрицы Tj наиболее подходящим методом. Полученный вектор Lj есть искомые собственные значения.
1.6 Последовательная сложность алгоритма
Все дальнейшие выкладки верны для наиболее быстрого последовательного варианта выполнения указываемых операций.
- Умножение квадратной матрицы порядка [math]n[/math] на вектор длины [math]n[/math] и наоборот требует по [math]n^2[/math] умножений и сложений.
- Перемножение векторов длины [math]n[/math] требует по [math]n[/math] умножений и сложений.
- Поэлементное сложение векторов длины [math]n[/math] требует [math]n[/math] сложений.
- Умножение вектора длины [math]n[/math] на число требует [math]n[/math] умножений.
- Нахождение квадратичной нормы вектора длины [math]n[/math] требует по [math]n[/math] умножений и сложений, а также одну операцию извлечения квадратного корня.
- Вычисление собственных значений матрицы порядка [math]k[/math] QR-алгоритмом требует [math]O(k^2)[/math] операций, вычисление также и собственных векторов требует примерно [math]6k ^ 3[/math], то есть [math]O(k^3)[/math] операций; при использовании метода «Разделяй-и-властвуй» для вычисления собственных значений и векторов аналогичной матрицы в среднем затрачивается [math]O(k^{2})[/math] операций.
Таким образом, для выполнения одной итерации метода Ланцоща требуется 1 операция вычисления квадратного корня, [math]n[/math] умножений, а также по [math]n ^ 2 + n + 2n + n[/math] сложений и умножений, а также суммарно [math]O(k ^ 3)[/math] операций, требуемых для поиска собственных значений и собственных векторов матрицы [math]T_{j}[/math]. То есть последовательная сложность алгоритма Ланцоша составляет [math]O(n ^ 2) + O(k ^ 3)[/math]. Это, очевидно, меньше, чем [math]O(n ^ 3)[/math] операций, требуемых в алгоритмах вычисления всех собственных значений произвольных симметричных матриц, в чём и заключается выгода метода Ланцоша.
1.7 Информационный граф
Информационный граф алгоритма, определяемый в [2] требует следующего замечания: подразумевается, что, например, умножение квадратной матрицы на вектор распараллеливается на [math]n[/math] потоков, где [math]n[/math] - порядок матрицы, а не на [math]n\log{n}[/math] потоков, чего можно было бы достичь, используя суммирование сдваиванием. Это допущение аналогично допущению из статьи Умножение матрицы на вектор и Это же допущение будет использовано и при оценке ресурса параллелизма алгоритма. Приводится граф без отображения входных и выходных данных.
1.8 Ресурс параллелизма алгоритма
Каждая итерация алгоритма Ланцоша выполняется строго последовательно. Однако внутри итерации возможно распараллеливание алгоритма.
- Умножение квадратной матрицы порядка [math]n[/math] на вектор длины [math]n[/math] требует последовательного выполнения [math]n[/math] ярусов умножений и сложений.
- Все остальные операции в рамках одной итерации, ведущие к вычислению матрицы [math]T_{j}[/math] выполняются строго последовательно. При этом вычисление значений векторов может быть выполнено за 1 ярус умножений / сложений, а вычисление чисел - за [math]\log{n}[/math] таких операций.
- Ресурс параллелизма алгоритма вычисления собственных значений зависит от используемого алгоритма и рассматривается в соответствующей статье.
Таким образом, при классификации по ширине ЯПФ алгоритм обладает линейной сложностью. В случае классификации по высоте ЯПФ алгоритм также имеет линейную сложность.
1.9 Входные и выходные данные алгоритма
Входные данные: квадратная матрица [math]A[/math] ( элементы [math]a_{ij}[/math] ), вектор [math]b[/math], число [math]k[/math]
Объём входных данных: [math]n^2 + n + 1[/math]
Выходные данные: вектор [math]L[/math], матрица [math]E[/math] ( элементы [math]e_{pq}[/math] )
Объём выходных данных: [math]k + kn[/math]
1.10 Свойства алгоритма
В случае неограниченности ресурсов соотношение последовательной и параллельной сложности алгоритма Ланцоща по сумме всех итераций представляет собой [math]k ^ 3 + kn ^ 2[/math] к [math]X + kn[/math], где [math]X[/math] - параллельная сложность алгоритма, используемого для вычисления собственных значений. В случае использования метода вычисления собственных значений с линейной сложностью при классисифкации по высоте ЯПФ это соотношение будет равно [math](k ^ 2 + n ^ 2) / (k + n)[/math].
Алгоритм Ланцоша не является полностью детерминированным в том смысле, что возможно выполнение числа итераций алгоритма, меньшего, чем заданное ( из-за того, что вычислены все ненулевые собственные значения матрицы [math]A[/math] ).
Важное свойство метода Ланцоша состоит в том, что первыми в матрице [math]T_{j}[/math] появляются собственные значения с максимальной величиной по модулю. Таким образом, метод особенно хорошо подходит для вычисления собственных значений матрицы [math]A[/math], находящихся на краях её спектра.
2 Программная реализация алгоритма
2.1 Особенности реализации последователього алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Для проверки масштабируемости алгоритма из реализации была исключена непосредственно задача поиска собственных значений матрицы [math]T_{j}[/math] на каждой итерации, так как это является предметом отдельных статей (QR-алгоритм, метод «Разделяй-и-властвуй», итерации Арнольди и.т.д.). Таким образом, проверяемая на масштабируемость задача сводится к итерационному вычислению коэффициентов матрицы [math]T_{j}[/math]. Матрица [math]A[/math] в рассматриваемой реализации хранится построчно и построчно же разделяется между процессорами. Также в целях исключения влияния недетерминированности алгоритма этот процесс изменен так, что при вычислении всех ненулевых собственных значений матрицы [math]A[/math] работа алгоритма продолжается над фиктивными данными, пока не будет завершено заранее заданное число итераций. При этом в реализации используется OpenMP для распараллеливания в пределах одного узла.
Набор изменяемых параметров запуска реализации алгоритма и границы значений параметров алгоритма:
- Число MPI процессов [1 : 200]
- Число OpenMP нитей [1 : 4] (Запуск с 1 OpenMP нитью производился только для 1 MPI процесса, ускорение замерялось относительно этого запуска).
- Размерность матрицы [400 : 4000]
Эффективность выполнения реализации алгоритма:
- Минимальная эффективность - <1% (неактуально, так как в исследовании достигнут предел масштабируемости алгоритма).
- Максимальная эффективность - ~30% (получается на 2-8 процессов при ~640000 ячеек матрицы на процесс; в случае 8 процессов это входная матрица 1600[math]{\cdot}[/math]1600).
Оценка масштабируемости:
- По числу процессов - при увеличении числа процессов эффективность уменьшается на всей области рассмотренных значений, но не слишком интенсивно. Это связано с небольшим объёмом пересылок данных, реализованных через достаточно эффективные операции редукции.
- По размеру задачи - при увеличении размера задачи эффективность вычислений вначале кратковременно возрастает, но затем начинает относительно равномерно убывать на всей рассматриваемой области.
- По двум направлениям - при рассмотрении увеличения, как вычислительной сложности, так и числа процессов по всей рассмотренной области значений наибольшая эффективность наблюдается при соответствии числа процессов и размера задачи. По этой "диагонали" эффективность вначале кратковременно возрастает, затем начинает убывать.
Интересно, что для разного числа процессоров положительные скачки производительности (зачастую - довольно заметные) имеют место на разных объёмах входных данных.
Далее приводятся графики ускорения программы для разных объёмов входных данных. Ускорение считалось относительно варианта с 1 MPI процессом без дальнейшнего распараллеливания посредством OpenMP. На оси абсцисс отмечается число MPI процессов, число OpenMP нитей во всех случаях равно четырём.
Используемая реализация алгоритма Реализация алгоритма гибридная, использует как MPI версии 2.2 (использование MPI_IN_PLACE), так и OpenMP. Указанный код компилировался с помощью компилятора IBM (mpixlcxx_r) с опцией -qsmp=omp (она включает в себя оптимизацию -O2). Реализация алгоритма полностью собственная и не использует встроенные библиотеки. Каждая "quadcore node" - это 1 узел системы BG/P, содержащий 4 микропроцессорных ядра IBM PowerPC 450 c суммарной пиковой производительностью 13.6 GFLOPS на узел и 2 GB общей памяти с пропускной способностью 13.6 GB/sec.
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
Подробный анализ алгоритмов нахождения собственных значений, включая алгоритм Ланцоша, был проведен в [3]. Существует несколько как последовательных, так и параллельных реализаций алгоритма Ланцоша, включенных в программные библиотеки для поиска собственных значений матриц. Свойства этих реализаций можно увидеть в таблицах ниже. Первая таблица — относительно недавние реализации, вторая - «музейные экспонаты». Следует уделить особенное внимание столбцам «тип метода» и «параллелизация». В первом из них значение только N соответствует описанному варианту без реортогонализации, F – полной реортогонализации, P – частичной реортогонализации, S – выборочной реортогонализации. Во втором значние «none» соответствует отсутствию параллельной реализации, M – параллельной реализации посредством MPI, O – параллельной реализации посредством OpenMP. Их приведение целесообразно, так как на практике алгоритм Ланцоша без переортогонализации неустойчив. Также указываются библиотеки, в которых реализованы более глубокие модификации метода Ланцоша, с указанием изменений в графе «тип метода».
| Название библиотеки | Язык программирования | Дата появления, версия библиотки | Тип метода | Параллелизация |
|---|---|---|---|---|
| BLKLAN | C/Matlab | 2003 | P | none |
| BLZPACK | F77 | 2000, 04/00 | P + S | M |
| IETL | C++ | 2006, 2.2 | N | none |
| SLEPc | C/F77 | 2009, 3.0.0 | All | M |
| TRLAN | F90 | 2006 | Dynamic thick-restart | M |
| PROPACK | F77/Matlab | 2005, 2.1 / 1.1 | P, finds SVD | O |
| IRBLEIGS | Matlab | 2000, 1.0 | Indefinitie symmetric | none |
| Название библиотеки | Язык программирования | Дата появления, версия библиотки | Тип метода | Параллелизация |
|---|---|---|---|---|
| ARPACK | F77 | 1995, 2.1 | Arnoldi iterations, impliicit restart | M |
| ARPACK++ | C++ | 1998, 1.1 | Arnoldi iterations, impliicit restart | none |
| LANCZOS | F77 | 1992 | N | none |
| LANZ | F77 | 1991, 1.0 | P | none |
| LASO | F77 | 1983, 2 | S | none |
| NAPACK | F77 | 1987 | N | none |
| QMRPACK | F77 | 1996 | Designed for nonsymmetrical matrices (lookahead) | none |
| SVDPACK | C/F77 | 1992 | P, finds SVD | none |
| Underwood | F77 | 1975 | F, block version | none |
3 Литература
<references \>
- ↑ Деммель Дж. Вычислительная линейная алгебра. Теория и приложения - М.: Мир, 2001. Ознакомиться
- ↑ Параллельные вычисления (Воеводин В.В., Воеводин Вл.В.) - Спб, изд-во "БХВ-Петербург", 2002 Ознакомиться
- ↑ A Survey of Software for Sparse Eigenvalue Problems - V. Hernandez, J. E. Roman, A. Tomas, V. Vidal (см.)
