Учacтник:Malikovmt/Алгоритм Ланцоша для арифметики с плавающей точкой с полной переортогонализацией: различия между версиями
[непроверенная версия] | [непроверенная версия] |
Malikovmt (обсуждение | вклад) |
ASA (обсуждение | вклад) |
||
(не показано 38 промежуточных версий 2 участников) | |||
Строка 1: | Строка 1: | ||
+ | {{Assignment|ASA}} | ||
+ | |||
{{algorithm | {{algorithm | ||
| name = Алгоритм Ланцоша с полной переортогонализацией | | name = Алгоритм Ланцоша с полной переортогонализацией | ||
| serial_complexity = <math>O(n^2k)</math> | | serial_complexity = <math>O(n^2k)</math> | ||
− | |||
− | |||
| input_data = <math>\frac{n(n + 1)}{2}</math> | | input_data = <math>\frac{n(n + 1)}{2}</math> | ||
| output_data = <math>k(n + 1)</math> | | output_data = <math>k(n + 1)</math> | ||
Строка 14: | Строка 14: | ||
=== Общее описание алгоритма === | === Общее описание алгоритма === | ||
− | '''Алгоритм Ланцоша''' - итерационный метод, используемый для вычисления части собственных значений и соответствующих им собственных векторов матрицы <math>A</math> размера <math>n*n</math>, изначально разработанный Корнелием Ланцошем. Преимуществами использования метода является относительно небольшое потребление памяти и вычислительных ресурсов, а также наличие параметра | + | '''Алгоритм Ланцоша''' - итерационный метод, используемый для вычисления части собственных значений и соответствующих им собственных векторов матрицы <math>A</math> размера <math>n*n</math>, изначально разработанный Корнелием Ланцошем. Преимуществами использования метода является относительно небольшое потребление памяти и вычислительных ресурсов, а также наличие параметра <math>k << n</math>, контролирующего количество итераций. Несмотря на то, что алгоритм является вычислительно эффективным, первоначально сформулированный метод был плохо применим из-за численной неустойчивости - метод хорошо работал на целочисленных значениях, однако в арифметике с плавающей точкой ошибки округления давали большую погрешность. В 1970 году Ojalvo и Newman показали, как сделать метод численно стабильным и применили его для расчета крупных инженерных сооружений, подверженных динамическим нагрузкам. Кроме того, они показали способ выбора начального приближения (с использованием ГПСЧ), а также эмпирический способ для выбора числа <math>k</math> (примерно в полтора раза больше искомого числа собственных векторов). В данный момент существует две основных модификации метода (с полной и выборочной переортогонализацией), а также большое количество модификаций, использующихся в различных технических областях. Алгоритм используется для больших <math>n</math>. |
=== Математическое описание алгоритма === | === Математическое описание алгоритма === | ||
+ | |||
+ | Первый этап алгоритма - использование метода Ланцоша для построения крыловского подпространства: <math> K_k(A,x) = span[x_1, Ax_1, A^2x_1, ..., A^{k-1}x_1] </math>. Входные данные алгоритма: квадратная симметричная матрица <math>A</math> размерности <math>n*n</math>, вектор начального приближения <math>b</math>, а так же число итераций <math>k</math>. Метод осуществляет поиск трехдиагональной симметричной матрицы <math>T_k=Q_k^TAQ_k</math>. | ||
+ | |||
+ | <math>T_k=\begin{bmatrix} | ||
+ | \alpha_1 & \beta_2 \\ | ||
+ | \beta_2 & \alpha_2 & \beta_3 &\\ | ||
+ | &. & . & .\\ | ||
+ | &&\beta_{k-1} & \alpha_{k-1} & \beta_k\\ | ||
+ | &&&\beta_k & \alpha_k | ||
+ | \end{bmatrix}</math> | ||
+ | |||
+ | Описание метода: | ||
+ | |||
+ | <math> | ||
+ | \begin{array}{l} | ||
+ | q_1 = b / \Vert b \Vert_2\\ | ||
+ | j = \overline{1, k}:\\ | ||
+ | \quad z_j = A q_j \\ | ||
+ | \quad \alpha_j = q_j^T z_j \\ | ||
+ | \quad z_j = z_j - \sum_{i=1}^j (z_j^T q_i) q_i\\ | ||
+ | \quad z_j = z_j - \sum_{i=1}^j (z_j^T q_i) q_i\\ | ||
+ | \quad \beta_j = \Vert z_j \Vert_2\\ | ||
+ | \quad q_{j+1} = z_j / \Vert z_j \Vert_2 = z_j/\beta_j | ||
+ | \end{array} | ||
+ | </math> | ||
+ | |||
+ | Следующий шаг алгоритма - процедура Рэлея-Ритца. Она зкалючается в интерпретации собственных значений матрицы <math> T_k=Q_k^TAQ_k</math>. | ||
+ | Ее собственные значения приближают собственные значения исходной матрицы. Пусть ''T<sub>k</sub>=V<math>\Lambda</math>V<sup>T</sup>'' - спектральное разложение матрицы ''T<sub>k</sub>'', тогда столбцы матрицы ''Q<sub>k</sub>V'' рассматриваются как приближения к соответствующим собственным векторам матрицы ''A'' и называются векторами Ритца. Числа и векторы Ритца являются оптимальными приближениями к собственным значениям и собственным векторам матрицы ''A''. | ||
+ | |||
+ | Поиск собственных значений матрицы ''T'' намного легче, чем для исходной матрицы, так как предполагается, что <math>k << n</math>, и матрица ''T'' - трехдиагональная. | ||
+ | |||
+ | Полная переортогонализация необходима для того, чтобы гарантировать, что каждый полученный вектор ''q<sub>j+1</sub>'' ортогонален уже имеющимся векторам ''q<sub>1..j</sub>''. Без этого процесса будут накапливаться существенные вычислительный ошибки. | ||
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
+ | |||
+ | Вычислительным ядро алгоритма состоит ииз двух основных частей: | ||
+ | * <math>Aq=( \sum\nolimits_{i=^n}a_{1i}q_i, \sum\nolimits_{i=2}^na_{2i}q_i, ..., \sum\nolimits_{i=1}^na_{ni}q_i)</math> - умножение симметричной матрицы <math>A</math> размерности <math>n*n</math> на вектор ''q'' размерности ''n''. | ||
+ | |||
+ | * <math>z=z-\sum\nolimits_{i=1}^{k}(z^Tq_i)q_i.</math> - процесс ортогонализации Грама-Шмидта. | ||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
+ | Основные операции алгоритма: | ||
+ | |||
+ | 1. Перемножение матрицы на вектор. | ||
+ | <math>b=Ax</math>. | ||
+ | |||
+ | 2. Двойная ортогонализация методом Грмма-Шмидта. | ||
+ | <math>z=z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i,</math> | ||
+ | <math>z=z-\alpha_jq_j-\beta_{j}q_{j-1}.</math> | ||
+ | |||
+ | 3. Вычисление обновленного базисного вектора. | ||
+ | <math>q_{i+1} = z/ \beta</math>. | ||
=== Схема реализации последовательного алгоритма === | === Схема реализации последовательного алгоритма === | ||
+ | |||
+ | В параграфе 1.2 приводится полная схема последовательного алгоритма. | ||
+ | |||
+ | Заполняем начальные значения алгоритма (''b'' - начальное преближение). | ||
+ | |||
+ | <math> | ||
+ | \begin{align} | ||
+ | & q_1=b/||b||,\\ | ||
+ | & \beta_1=0,\\ | ||
+ | &q_0=0, \\ | ||
+ | \end{align} | ||
+ | </math> | ||
+ | |||
+ | Для всех <math>j=1..k</math>: | ||
+ | |||
+ | *1. Вычисляется ''j''-й диагональный элемент матрицы <math>T_k</math>: <math>z=Aq_j; \alpha_j=q_j^Tz;</math> | ||
+ | |||
+ | *2. Проводится полная переортогонализация Грамма-Шмидта: <math>z =z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i;</math> | ||
+ | |||
+ | *3. Вычисляются значения <math>\beta_{j+1}</math> матрицы <math>T_k</math>: <math>\beta_{j+1}=||z||;</math> | ||
+ | |||
+ | *4. Если <math>\beta_{j+1}=0</math>, то алгоритм завершается; | ||
+ | |||
+ | *5. Сохраняем значения для следующей итерации <math>q_{j+1}=z/\beta_{j+1}.</math> | ||
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === | ||
+ | |||
+ | *1. Основная часть операций в алгоритме Ланцоша производится во время умножения матрицы <math>A</math> размерности <math>n*n</math> на вектор <math>q</math> размерности <math>n</math> - вычислительная сложность: <math>n^2</math> умножений и <math>n^2-n</math> сложений. Остальные операции основного цикла производят меньше <math>n^2</math> операций сложения или умножения. Так как умножение матрицы на вектор производится <math>k</math> раз, то сложность этой части алгоритма - <math>O(kn^2)</math> | ||
+ | *2. Процесс ортогонализации Грама-Шмидта - вычислительная сложность: <math>k^2n+k(n+2)</math> умножений и <math>k^2n + k(n + 1) + 2</math> сложений. Производится в цикле <math>k</math> раз. Сложность - <math>O(nk^2)</math> | ||
+ | *3. Процесс разложения матрицы <math>T</math> размерности <math>k*k</math>. Сложность - <math>O(k^2)</math> | ||
+ | |||
+ | Так как число итераций много меньше размерности матрицы <math>A</math>, <math>k << n</math>, то общая сложность алгоритма сокращается до <math>O(kn^2)</math>. | ||
=== Информационный граф === | === Информационный граф === | ||
+ | |||
+ | [[Файл:Lanczos.png|center|300px|'''Рис. 1.''' Граф ''j''-й итерации.]] | ||
+ | |||
+ | Черные маленькие стрелочки - куда идут данные. | ||
+ | |||
+ | Большие стрелочки - обозначение прохода по массиву (изменение индекса) | ||
+ | |||
+ | В больших кружках преобразования (операции над данными или вычисления) | ||
+ | |||
+ | В маленьких кружках запись в переменную. | ||
=== Ресурс параллелизма алгоритма === | === Ресурс параллелизма алгоритма === | ||
+ | |||
+ | Алгоритм Ланцоша - итерационный, итерации должны выполняться в строгой последовательности, и нет возможности их параллелизовать. Внутри одной итерации алгоритма ресурсами параллелизма могут быть: | ||
+ | *1. умножение матрицы размерности <math>n * n</math> на вектор длины <math>n</math> требует последовательного выполнения <math>n</math> ярусов умножений и сложений; | ||
+ | *2. вычисление <math>\alpha_j</math> требует <math>n</math> ярусов сложений с 1 операцией умножения в каждом; | ||
+ | *3. переортагонализация требует вычисления <math>j</math> ярусов сложений с <math>n</math> операциями умножения в каждом, <math>n</math> ярусов сложений с <math>j</math> операциями умножения в каждом; | ||
+ | *4. вычисление нормы вектора длины <math>n</math> требует <math>n</math> ярусов сложений с <math>1</math> операцией умножения в каждом. | ||
=== Входные и выходные данные алгоритма === | === Входные и выходные данные алгоритма === | ||
+ | '''Входные данные''': <math>A\in\mathbb{R}^{n\times n}</math> - симметричная матрица, т. е. <math>a_{ij}= a_{ji}, i, j = 1, \ldots, n</math>. Объем данных: <math>\frac{n (n + 1)}{2}</math>. | ||
+ | |||
+ | '''Выходные данные''': <math>\Lambda</math> - вектор собственных значений матрицы и соответствующие им собственные вектора <math>v_i,i=1,\ldots, k </math>. Объем данных: <math> k(n+1) </math>. | ||
=== Свойства алгоритма === | === Свойства алгоритма === | ||
+ | * В классическом алгоритме Ланцоша возникает большая погрешность при округлении чисел с плавающей точкой. Выбранный вариант с полной переортогонализацией устраняет этот недостаток, однако является более ресурсоемким. На практике наиболее популярен вариант с частичной переортогонализацией. | ||
+ | * Вычислительная мощность алгоритма (отношение числа операций к суммарному объему данных) оценивается как <math>\approx 2k</math>. | ||
+ | * Преимуществом алгоритма является то, что он начинает поиск собственных значений матрицы начиная с максимального в абсолютном смысле значения. | ||
+ | * Нет необходимости хранения исходной матрицы на каждом вычислительном ядре, так как метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать регулярность структуры матрицы. | ||
+ | * Существует модифицированный блочный вариант алгоритма, применяемый в случае кратных собственных значений. | ||
+ | |||
+ | == Программная реализация алгоритма == | ||
+ | |||
+ | === Особенности реализации последовательного алгоритма === | ||
=== Локальность данных и вычислений === | === Локальность данных и вычислений === | ||
+ | === Возможные способы и особенности параллельной реализации алгоритма === | ||
+ | === Масштабируемость алгоритма и его реализации === | ||
− | + | Проверка масштабируемости алгоритма проходила на [//hpc.cmc.msu.ru/ IBM Blue Gene/P ВМК МГУ.] | |
− | + | Используемые компиляторы: intel/15.0.090 и OpenMPI/1.8.4-icc. | |
− | + | Для проведения расчетов и получения полноценной картины поведения алгоритма в зависимости от входных данных и числа процессоров, программа была запущена на следующих параметрах: | |
+ | * размеры входной матрицы: [5000:30000] c шагом 5000; | ||
+ | * число процессоров: от 1, 2, 4 .. 64. | ||
− | + | Сборка осуществлялась с параметрами: | |
+ | *openmpi/1.5.5-icc | ||
+ | *intel/13.1.0 | ||
− | + | [[Файл:Lancos prof.jpg|центр|мини|800x800пкс|Рис. 2. Зависимость времени выполнения программы от входных данных]] | |
− | + | По графику видно, что время выполнения уменьшается с увеличением количества процессоров примерно до 16, дальше начинает медленно увеличиваться вплоть до 64 и далее. Это можно объяснить тем, что с увеличением количества вычислительных ядер растет количество передаваемых данных, и дальнейшее увеличение их числа только ухудшает положение. | |
− | + | [http://pastebin.com/sWQ9ihcW код программы] | |
=== Динамические характеристики и эффективность реализации алгоритма === | === Динамические характеристики и эффективность реализации алгоритма === | ||
+ | === Выводы для классов архитектур === | ||
+ | === Существующие реализации алгоритма === | ||
+ | |||
+ | The IETL Project http://www.comp-phys.org/software/ietl/ C++ | ||
+ | |||
+ | NAG Library http://www.nag.com/content/nag-library C, C++, Fortran, C#, MATLAB, R | ||
+ | |||
+ | ARPACK https://people.sc.fsu.edu/~jburkardt/m_src/arpack/arpack.html MATLAB | ||
+ | |||
+ | GrapLab https://turi.com/products/create/open_source.html C++ | ||
− | + | LANSO/PLANSO http://web.cs.ucdavis.edu/~bai/ET/lanczos_methods/overview_PLANSO.html Fortran (уже распараллелена) | |
+ | |||
+ | Julia Math https://github.com/JuliaMath/IterativeSolvers.jl Julia | ||
− | + | SciPy https://scipy.org/ Python | |
== Литература == | == Литература == | ||
+ | |||
+ | Дж. Деммель «Вычислительная линейная алгебра» (стр. 391) | ||
<references \> | <references \> | ||
− | |||
− | |||
− | |||
#[[Категория:Разложения матриц]] | #[[Категория:Разложения матриц]] |
Текущая версия на 13:15, 4 октября 2017
Эта работа прошла предварительную проверку Дата последней правки страницы: 04.10.2017 Данная работа соответствует формальным критериям. Проверено ASA. |
Алгоритм Ланцоша с полной переортогонализацией | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(n^2k)[/math] |
Объём входных данных | [math]\frac{n(n + 1)}{2}[/math] |
Объём выходных данных | [math]k(n + 1)[/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]A[/math] размера [math]n*n[/math], изначально разработанный Корнелием Ланцошем. Преимуществами использования метода является относительно небольшое потребление памяти и вычислительных ресурсов, а также наличие параметра [math]k \lt \lt n[/math], контролирующего количество итераций. Несмотря на то, что алгоритм является вычислительно эффективным, первоначально сформулированный метод был плохо применим из-за численной неустойчивости - метод хорошо работал на целочисленных значениях, однако в арифметике с плавающей точкой ошибки округления давали большую погрешность. В 1970 году Ojalvo и Newman показали, как сделать метод численно стабильным и применили его для расчета крупных инженерных сооружений, подверженных динамическим нагрузкам. Кроме того, они показали способ выбора начального приближения (с использованием ГПСЧ), а также эмпирический способ для выбора числа [math]k[/math] (примерно в полтора раза больше искомого числа собственных векторов). В данный момент существует две основных модификации метода (с полной и выборочной переортогонализацией), а также большое количество модификаций, использующихся в различных технических областях. Алгоритм используется для больших [math]n[/math].
1.2 Математическое описание алгоритма
Первый этап алгоритма - использование метода Ланцоша для построения крыловского подпространства: [math] K_k(A,x) = span[x_1, Ax_1, A^2x_1, ..., A^{k-1}x_1] [/math]. Входные данные алгоритма: квадратная симметричная матрица [math]A[/math] размерности [math]n*n[/math], вектор начального приближения [math]b[/math], а так же число итераций [math]k[/math]. Метод осуществляет поиск трехдиагональной симметричной матрицы [math]T_k=Q_k^TAQ_k[/math].
[math]T_k=\begin{bmatrix} \alpha_1 & \beta_2 \\ \beta_2 & \alpha_2 & \beta_3 &\\ &. & . & .\\ &&\beta_{k-1} & \alpha_{k-1} & \beta_k\\ &&&\beta_k & \alpha_k \end{bmatrix}[/math]
Описание метода:
[math] \begin{array}{l} q_1 = b / \Vert b \Vert_2\\ j = \overline{1, k}:\\ \quad z_j = A q_j \\ \quad \alpha_j = q_j^T z_j \\ \quad z_j = z_j - \sum_{i=1}^j (z_j^T q_i) q_i\\ \quad z_j = z_j - \sum_{i=1}^j (z_j^T q_i) q_i\\ \quad \beta_j = \Vert z_j \Vert_2\\ \quad q_{j+1} = z_j / \Vert z_j \Vert_2 = z_j/\beta_j \end{array} [/math]
Следующий шаг алгоритма - процедура Рэлея-Ритца. Она зкалючается в интерпретации собственных значений матрицы [math] T_k=Q_k^TAQ_k[/math]. Ее собственные значения приближают собственные значения исходной матрицы. Пусть Tk=V[math]\Lambda[/math]VT - спектральное разложение матрицы Tk, тогда столбцы матрицы QkV рассматриваются как приближения к соответствующим собственным векторам матрицы A и называются векторами Ритца. Числа и векторы Ритца являются оптимальными приближениями к собственным значениям и собственным векторам матрицы A.
Поиск собственных значений матрицы T намного легче, чем для исходной матрицы, так как предполагается, что [math]k \lt \lt n[/math], и матрица T - трехдиагональная.
Полная переортогонализация необходима для того, чтобы гарантировать, что каждый полученный вектор qj+1 ортогонален уже имеющимся векторам q1..j. Без этого процесса будут накапливаться существенные вычислительный ошибки.
1.3 Вычислительное ядро алгоритма
Вычислительным ядро алгоритма состоит ииз двух основных частей:
- [math]Aq=( \sum\nolimits_{i=^n}a_{1i}q_i, \sum\nolimits_{i=2}^na_{2i}q_i, ..., \sum\nolimits_{i=1}^na_{ni}q_i)[/math] - умножение симметричной матрицы [math]A[/math] размерности [math]n*n[/math] на вектор q размерности n.
- [math]z=z-\sum\nolimits_{i=1}^{k}(z^Tq_i)q_i.[/math] - процесс ортогонализации Грама-Шмидта.
1.4 Макроструктура алгоритма
Основные операции алгоритма:
1. Перемножение матрицы на вектор. [math]b=Ax[/math].
2. Двойная ортогонализация методом Грмма-Шмидта. [math]z=z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i,[/math] [math]z=z-\alpha_jq_j-\beta_{j}q_{j-1}.[/math]
3. Вычисление обновленного базисного вектора. [math]q_{i+1} = z/ \beta[/math].
1.5 Схема реализации последовательного алгоритма
В параграфе 1.2 приводится полная схема последовательного алгоритма.
Заполняем начальные значения алгоритма (b - начальное преближение).
[math] \begin{align} & q_1=b/||b||,\\ & \beta_1=0,\\ &q_0=0, \\ \end{align} [/math]
Для всех [math]j=1..k[/math]:
- 1. Вычисляется j-й диагональный элемент матрицы [math]T_k[/math]: [math]z=Aq_j; \alpha_j=q_j^Tz;[/math]
- 2. Проводится полная переортогонализация Грамма-Шмидта: [math]z =z-\sum\nolimits_{i=1}^{j-1}(z^Tq_i)q_i;[/math]
- 3. Вычисляются значения [math]\beta_{j+1}[/math] матрицы [math]T_k[/math]: [math]\beta_{j+1}=||z||;[/math]
- 4. Если [math]\beta_{j+1}=0[/math], то алгоритм завершается;
- 5. Сохраняем значения для следующей итерации [math]q_{j+1}=z/\beta_{j+1}.[/math]
1.6 Последовательная сложность алгоритма
- 1. Основная часть операций в алгоритме Ланцоша производится во время умножения матрицы [math]A[/math] размерности [math]n*n[/math] на вектор [math]q[/math] размерности [math]n[/math] - вычислительная сложность: [math]n^2[/math] умножений и [math]n^2-n[/math] сложений. Остальные операции основного цикла производят меньше [math]n^2[/math] операций сложения или умножения. Так как умножение матрицы на вектор производится [math]k[/math] раз, то сложность этой части алгоритма - [math]O(kn^2)[/math]
- 2. Процесс ортогонализации Грама-Шмидта - вычислительная сложность: [math]k^2n+k(n+2)[/math] умножений и [math]k^2n + k(n + 1) + 2[/math] сложений. Производится в цикле [math]k[/math] раз. Сложность - [math]O(nk^2)[/math]
- 3. Процесс разложения матрицы [math]T[/math] размерности [math]k*k[/math]. Сложность - [math]O(k^2)[/math]
Так как число итераций много меньше размерности матрицы [math]A[/math], [math]k \lt \lt n[/math], то общая сложность алгоритма сокращается до [math]O(kn^2)[/math].
1.7 Информационный граф
Черные маленькие стрелочки - куда идут данные.
Большие стрелочки - обозначение прохода по массиву (изменение индекса)
В больших кружках преобразования (операции над данными или вычисления)
В маленьких кружках запись в переменную.
1.8 Ресурс параллелизма алгоритма
Алгоритм Ланцоша - итерационный, итерации должны выполняться в строгой последовательности, и нет возможности их параллелизовать. Внутри одной итерации алгоритма ресурсами параллелизма могут быть:
- 1. умножение матрицы размерности [math]n * n[/math] на вектор длины [math]n[/math] требует последовательного выполнения [math]n[/math] ярусов умножений и сложений;
- 2. вычисление [math]\alpha_j[/math] требует [math]n[/math] ярусов сложений с 1 операцией умножения в каждом;
- 3. переортагонализация требует вычисления [math]j[/math] ярусов сложений с [math]n[/math] операциями умножения в каждом, [math]n[/math] ярусов сложений с [math]j[/math] операциями умножения в каждом;
- 4. вычисление нормы вектора длины [math]n[/math] требует [math]n[/math] ярусов сложений с [math]1[/math] операцией умножения в каждом.
1.9 Входные и выходные данные алгоритма
Входные данные: [math]A\in\mathbb{R}^{n\times n}[/math] - симметричная матрица, т. е. [math]a_{ij}= a_{ji}, i, j = 1, \ldots, n[/math]. Объем данных: [math]\frac{n (n + 1)}{2}[/math].
Выходные данные: [math]\Lambda[/math] - вектор собственных значений матрицы и соответствующие им собственные вектора [math]v_i,i=1,\ldots, k [/math]. Объем данных: [math] k(n+1) [/math].
1.10 Свойства алгоритма
- В классическом алгоритме Ланцоша возникает большая погрешность при округлении чисел с плавающей точкой. Выбранный вариант с полной переортогонализацией устраняет этот недостаток, однако является более ресурсоемким. На практике наиболее популярен вариант с частичной переортогонализацией.
- Вычислительная мощность алгоритма (отношение числа операций к суммарному объему данных) оценивается как [math]\approx 2k[/math].
- Преимуществом алгоритма является то, что он начинает поиск собственных значений матрицы начиная с максимального в абсолютном смысле значения.
- Нет необходимости хранения исходной матрицы на каждом вычислительном ядре, так как метод использует исходную матрицу только в операциях умножения матрицы на вектор, что позволяет эффективно использовать регулярность структуры матрицы.
- Существует модифицированный блочный вариант алгоритма, применяемый в случае кратных собственных значений.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Проверка масштабируемости алгоритма проходила на IBM Blue Gene/P ВМК МГУ.
Используемые компиляторы: intel/15.0.090 и OpenMPI/1.8.4-icc.
Для проведения расчетов и получения полноценной картины поведения алгоритма в зависимости от входных данных и числа процессоров, программа была запущена на следующих параметрах:
- размеры входной матрицы: [5000:30000] c шагом 5000;
- число процессоров: от 1, 2, 4 .. 64.
Сборка осуществлялась с параметрами:
- openmpi/1.5.5-icc
- intel/13.1.0
По графику видно, что время выполнения уменьшается с увеличением количества процессоров примерно до 16, дальше начинает медленно увеличиваться вплоть до 64 и далее. Это можно объяснить тем, что с увеличением количества вычислительных ядер растет количество передаваемых данных, и дальнейшее увеличение их числа только ухудшает положение.
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации алгоритма
The IETL Project http://www.comp-phys.org/software/ietl/ C++
NAG Library http://www.nag.com/content/nag-library C, C++, Fortran, C#, MATLAB, R
ARPACK https://people.sc.fsu.edu/~jburkardt/m_src/arpack/arpack.html MATLAB
GrapLab https://turi.com/products/create/open_source.html C++
LANSO/PLANSO http://web.cs.ucdavis.edu/~bai/ET/lanczos_methods/overview_PLANSO.html Fortran (уже распараллелена)
Julia Math https://github.com/JuliaMath/IterativeSolvers.jl Julia
SciPy https://scipy.org/ Python
3 Литература
Дж. Деммель «Вычислительная линейная алгебра» (стр. 391)
<references \>