Участник:A.Freeman/Алгоритм Ланцоша для точной арифметики (без переортогонализации): различия между версиями
A.Freeman (обсуждение | вклад) |
ASA (обсуждение | вклад) |
||
(не показаны 52 промежуточные версии 4 участников) | |||
Строка 1: | Строка 1: | ||
+ | {{Assignment|ASA|Konshin}} | ||
+ | |||
{{algorithm | {{algorithm | ||
| name = Алгоритм Ланцоша без переортогонализации | | name = Алгоритм Ланцоша без переортогонализации | ||
Строка 7: | Строка 9: | ||
| output_data = <math>k(n+1)</math> | | output_data = <math>k(n+1)</math> | ||
}} | }} | ||
− | + | Основные авторы описания: [[Участник:KAKTUZ|Заспа А.Ю.]] ([[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Вычислительное ядро алгоритма|1.3]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Информационный граф|1.7]]), [[Участник:A.Freeman|Фролов А.А.]] ([[#Последовательная сложность алгоритма|1.6]], [[#Информационный граф|1.7]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]], [[#Масштабируемость алгоритма и его реализации|2.4]], [[#Существующие реализации|2.7]]) | |
− | Основные авторы описания: [[Участник:KAKTUZ|Заспа А.Ю.]] ([[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Вычислительное ядро алгоритма|1.3]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Информационный граф|1.7]]), [[Участник:A.Freeman|Фролов А.А.]] ([[#Последовательная сложность алгоритма|1.6]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Входные и выходные данные алгоритма|1.9]], [[#Свойства алгоритма|1.10]], [[#Существующие реализации|2.7]]) | ||
== Свойства и структура алгоритма == | == Свойства и структура алгоритма == | ||
Строка 38: | Строка 39: | ||
Вычисляемые данные: собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. | Вычисляемые данные: собственные вектора матрицы <math>T_k</math> являющиеся столбцами матрицы <math>Q_k V</math>, и матрица собственных значений <math>\Lambda</math>, где <math>V, \Lambda</math> из спектрального разложения <math>T_k = V\Lambda V^T</math>. | ||
+ | До начала первой итерации нормируется вектор начального приближения <math>q_1 = \frac{b}{\|b\|_2}</math>. Затем задаются константы <math>\beta_0 = 0,\; q_0 = 0 </math>. | ||
− | + | После выполнения данных действий начинаются итерации алгоритма. Которых будет не больше чем <math>k</math>. В ходе итераций постепенно формируется матрица <math>T_k</math> из <math>\alpha_i</math>-x и <math>\beta_i</math>-х. | |
− | <math> | + | На <math>i</math>-й итерации происходят следующие вычисления. Вначале вычисляется произведение матрицы на вектор <math>z = Aq_i</math>. Затем считается скалярное произведение <math>\alpha_i = q^T_i z</math> и значение сохраняется как элемент формируемой матрицы. После этого необходимо произвести вычисление линейной комбинации векторов <math>z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}</math>. Норму полученного вектора запишем как элемент матрицы <math>\beta_i = \|z\|_2</math>. Если <math>\beta_i</math> оказалась равной нулю, то больше никаких итераций в алгоритме не производят, сразу переходят к вычислению собственных значений и собственных векторов. Иначе, в конце итерации считают <math>q_{i+1} = \frac{z}{\beta_i}</math> и переходят к следующей итерации. |
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | |||
− | </math> | ||
− | + | После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы <math>T_k</math>.<ref>Деммель Д. Вычислительная линейная алгебра</ref> | |
=== Вычислительное ядро алгоритма === | === Вычислительное ядро алгоритма === | ||
Строка 64: | Строка 54: | ||
=== Макроструктура алгоритма === | === Макроструктура алгоритма === | ||
+ | Для удобства понимания алгоритма можно выделить следующие макрооперации: | ||
+ | * умножение матрицы на вектор. Состоит из операций умножения вектора на число и сложения векторов. | ||
+ | * скалярное произведение векторов. | ||
+ | * линейная комбинация векторов. Состоит из умножений вектора на число и сложений векторов. | ||
+ | * получение нормы вектора. Состоит из скалярного произведения векторов, вычисления квадратного корня. | ||
+ | * деление вектора на число. | ||
+ | * вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы. | ||
+ | |||
+ | На каждой итерации будет вычисляться умножение матрицы <math>A</math> размера <math>n \times n</math> на вектор <math>q_i</math>. Эта операция и будет самой ресурсозатратной на каждой итерации. Результаты этой макрооперации будут использованы затем в скалярном произведении векторов и в линейной комбинации из 3 векторов. После этого результат линейной комбинации векторов необходимо нормировать. Для этого можно разбить эту процедуру на вычисление нормы вектора и деление вектора на получившееся число. Кроме того, после завершения итераций необходимо вычислить собственные вектора и собственные значения полученной матрицы. Данная операция не конкретизируется алгоритмом, и может быть выполнена любым способом. Поэтому удобнее её представить в виде макрооперации. В разделе [[#Информационный граф|1.7]] можно наглядно увидеть как передаются данные между этими макрооперациями. | ||
+ | |||
− | |||
− | |||
− | |||
− | |||
− | |||
=== Схема реализации последовательного алгоритма === | === Схема реализации последовательного алгоритма === | ||
Строка 75: | Строка 70: | ||
Последовательность исполнения метода следующая: | Последовательность исполнения метода следующая: | ||
− | <math> | + | <math>1.\, \beta_0 = 0,\; q_0 = 0</math> #Инициализируются константы |
− | + | ||
− | 1. \, \beta_0 = | + | <math>2.\, \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}</math> #Вычисляем норму вектора начального приближения. |
− | \|b\|_2 = | + | |
− | + | <math>3.\, q_{1_{j}} = \frac{b_{j}}{\|b\|_2}, \; j = 1,\, \dots\,, n</math> #Нормализуем вектор начального приближения. | |
− | + | ||
− | + | Начинаются итерации цикла, которых не больше чем <math>k</math>. На <math>i</math>-й итерации производятся следующие вычисления: | |
− | + | ||
− | + | <math>i.1\, z_j = \sum\limits_{m=1}^{n} a_{jm} q_{i_m}, \; j = 1,\,\dots\,, n</math> #Считаем результат применения линейного оператора <math>A</math> к вектору <math>q_i</math>. | |
− | + | ||
− | + | <math>i.2\, \alpha_i = \sum\limits_{j=1}^{n}q_{i_j} z_j</math> #Получаем результат скалярного произведения векторов <math>q_i</math> и <math>z</math>. | |
− | + | ||
− | + | <math>i.3\, z_j = z_j - \alpha_i q_{i_j} - \beta_{i-1}q_{i-1_j}, \, j = 1,\,\dots\,, n</math> #Вычисляем линейную комбинацию векторов. | |
− | + | ||
− | </math> | + | <math>i.4\, \beta_i = \|z\|_2 = \sqrt{\sum\limits_{j=1}^{n} z_j^2}</math> #Считаем норму вектора <math>z</math>. |
+ | |||
+ | <math>i.5</math> Проверка равенства <math>\beta_i == 0</math> # Если норма оказалась равной нулю, то завершаем итерации и переходим к вычислению собственных векторов и собственных значений полученной матрицы. В обратном случае, продолжаем выполнения итераций. | ||
+ | |||
+ | <math>i.6\, q_{i+1_j} = \frac{z_j}{\beta_i}, \; j = 1,\, \dots \,, n</math> #Нормируем вектор <math>z</math>. | ||
+ | |||
+ | <math>i.7\,</math> Если выполнили <math>k</math> итераций, то завершаем выполнение итераций, переходим к следующему шагу. Иначе начинаем последующую итерацию цикла. | ||
+ | |||
+ | <math>4.</math> Вычисляем собственные значения и собственные вектора полученной матрицы <math>T_k</math>. | ||
=== Последовательная сложность алгоритма === | === Последовательная сложность алгоритма === | ||
Строка 97: | Строка 100: | ||
Инициализация состоит из следующих операций: | Инициализация состоит из следующих операций: | ||
* вычисление нормы: | * вычисление нормы: | ||
− | ** <math>n | + | ** <math>n</math> операций умножения |
** <math>n-1</math> операций сложений | ** <math>n-1</math> операций сложений | ||
** извлечение квадратного корня | ** извлечение квадратного корня | ||
Строка 118: | Строка 121: | ||
** <math>n</math> операций деления | ** <math>n</math> операций деления | ||
− | Итого (без учета решения трехдиагональной матрицы): | + | Итого (без учета решения задачи собственных значений трехдиагональной матрицы): |
* <math> k(n^2 + 4n) + n </math> умножений, | * <math> k(n^2 + 4n) + n </math> умножений, | ||
* <math> k(n^2+3n-2) + n-1</math> сложений/вычитаний, | * <math> k(n^2+3n-2) + n-1</math> сложений/вычитаний, | ||
Строка 130: | Строка 133: | ||
=== Информационный граф === | === Информационный граф === | ||
− | [[Файл:Ланцош без переортогонализации 1.png|600px|thumb|center|Рисунок 1. Граф алгоритма без отображения входных и выходных данных.<br/> | + | {|align="center" |
− | <math>\|x\|</math> — макрооперация взятия нормы, <math>/</math> — операция деления, | + | |-valign="top" |
− | <math>\mathsf{F}</math> — вычисление линейной комбинации векторов, | + | |[[Файл:Ланцош без переортогонализации 1.png|600px|thumb|center|Рисунок 1. Граф алгоритма без отображения входных и выходных данных.<br/> |
− | [[Файл:Ланцош без переортогонализации 2.png|600px|thumb|center|Рисунок 2. Граф алгоритма с отображением входных и выходных данных.<br/> | + | <math>\|x\|</math> — макрооперация взятия нормы, <br/> |
− | <math>\|x\|</math> — макрооперация взятия нормы, <math>/</math> — операция деления, <math>\cdot </math> — операция умножения,<br/> | + | <math>/</math> — операция деления, <br/> |
− | <math>\mathsf{F}</math> — вычисление линейной комбинации векторов, <math>\mathsf{Ax}</math> — линейный | + | <math>\cdot </math> — операция скалярного умножения векторов, <br/> |
− | <math>\mathsf{In}</math> — входные данные, <math>\mathsf{Out}</math> — выходные результат.]] | + | <math>\mathsf{F}</math> — вычисление линейной комбинации векторов, <br/> |
+ | <math>\mathsf{Ax}</math> — линейный оператор, <br/> | ||
+ | <math>\nabla</math> — условие преждевременного выхода.]] | ||
+ | |[[Файл:Ланцош без переортогонализации 2.png|600px|thumb|center|Рисунок 2. Граф алгоритма с отображением входных и выходных данных.<br/> | ||
+ | <math>\|x\|</math> — макрооперация взятия нормы, <br/> | ||
+ | <math>/</math> — операция деления, <br/> | ||
+ | <math>\cdot </math> — операция скалярного умножения векторов,<br/> | ||
+ | <math>\mathsf{F}</math> — вычисление линейной комбинации векторов, <br/> | ||
+ | <math>\mathsf{Ax}</math> — линейный оператор, <br/> | ||
+ | <math>\nabla</math> — условие преждевременного выхода,<br/> | ||
+ | <math>\mathsf{In}</math> — входные данные, <br/> | ||
+ | <math>\mathsf{Out}</math> — выходные результат.]] | ||
+ | |} | ||
+ | [[Файл:Линейный_оператор.png|600px|thumb|center|Рисунок 3. Граф макрооперации "умножение матрицы на вектор" без отображения входных и выходных данных.<br/> | ||
+ | <math>*</math> — операция умножения двух коэффициентов, <br/> | ||
+ | <math>+</math> — операция сложения.]] | ||
=== Ресурс параллелизма алгоритма === | === Ресурс параллелизма алгоритма === | ||
Строка 194: | Строка 212: | ||
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является субквадратичным: <math>\frac{kn^2}{k \log{n}}</math>. | Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является субквадратичным: <math>\frac{kn^2}{k \log{n}}</math>. | ||
+ | |||
+ | Вычислительная мощность алгоритма нахождения коэффициентов <math>\alpha_i,\beta_i</math> без экономии памяти состовляет <math>\frac{k(2n^2+8n-1)+3n}{n^2+2k}</math>, при <math>k \ll n</math>, вычислительная мощность приблизительно равна <math>2k</math>. | ||
+ | |||
+ | В случае хранения половины исходной матрицы получаем соотношношение <math>\frac{k(2n^2+8n-1)+3n}{n(n+1)/2+2k}</math>, что приблизительно равно <math>4k</math>. | ||
Алгоритм Ланцоша без переортогонализации не является детерминированным (возможно выполнение меньшего числа итераций алгоритма), в случае если все собственные значения уже вычислены. | Алгоритм Ланцоша без переортогонализации не является детерминированным (возможно выполнение меньшего числа итераций алгоритма), в случае если все собственные значения уже вычислены. | ||
Строка 208: | Строка 230: | ||
=== Масштабируемость алгоритма и его реализации === | === Масштабируемость алгоритма и его реализации === | ||
+ | Проведём исследование масштабируемости параллельной реализации алгоритма Ланцоша без переортогонализации согласно [[Scalability methodology|методике]]. Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета]. | ||
+ | |||
+ | ==== Реализация на языке C ==== | ||
+ | [https://github.com/xrstalker/lanczos Исследованная параллельная реализация на языке C] | ||
+ | |||
+ | ==== Масштабируемость реализации алгоритма ==== | ||
+ | |||
+ | Сборка осуществлялась со следующими параметрами: | ||
+ | |||
+ | * gcc-5.2.0 | ||
+ | * openmpi-1.8.4 | ||
+ | * аргументы компилятора: -std=c11 -Ofast | ||
+ | |||
+ | Набор и границы значений изменяемых параметров реализации алгоритма: | ||
+ | |||
+ | * число процессоров [32 : 256] с шагом 16; | ||
+ | * размер матрицы [5000 : 200000] с шагом 5000. | ||
+ | |||
+ | В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма: | ||
+ | |||
+ | * минимальная эффективность реализации 1.92%; | ||
+ | * максимальная эффективность реализации 59.47%. | ||
+ | |||
+ | График полученного распределения эффективности: | ||
+ | |||
+ | [[file:Lanczos Without Orthogonalization Eff2Dist.png|thumb|center|600px|Рисунок 4. Параллельная реализация алгоритма Ланцоша. Распределение производительности.]] | ||
+ | |||
+ | На следующих рисунках приведены графики производительности и эффективности данной реализации в зависимости от изменяемых параметров запуска. | ||
+ | {|align="center" | ||
+ | |-valign="top" | ||
+ | |[[file:Lanczos Without Orthogonalization Perf2.png|thumb|center|600px|Рисунок 5. Параллельная реализация алгоритма Ланцоша. Изменение производительности в зависимости от числа процессоров и размера матрицы.]] | ||
+ | |[[file:Lanczos Without Orthogonalization Eff2.png|thumb|center|600px|Рисунок 6. Параллельная реализация алгоритма Ланцоша. Изменение эффективности в зависимости от числа процессоров и размера матрицы.]] | ||
+ | |} | ||
+ | |||
+ | Средняя эффективность для матриц с размером больше 25000 составила 52.7%. | ||
=== Динамические характеристики и эффективность реализации алгоритма === | === Динамические характеристики и эффективность реализации алгоритма === |
Текущая версия на 16:58, 19 декабря 2016
Эта работа успешно выполнена Преподавателю: в основное пространство, в подстраницу Данное задание было проверено и зачтено. Проверено Konshin и ASA. |
Алгоритм Ланцоша без переортогонализации | |
Последовательный алгоритм | |
Последовательная сложность | [math]O(kn^2)[/math] |
Объём входных данных | [math]\frac{n (n + 1)}{2}+1[/math] |
Объём выходных данных | [math]k(n+1)[/math] |
Параллельный алгоритм | |
Высота ярусно-параллельной формы | [math]O(k \log{n})[/math] |
Ширина ярусно-параллельной формы | [math]O(n^2)[/math] |
Основные авторы описания: Заспа А.Ю. (1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7), Фролов А.А. (1.6, 1.7, 1.8, 1.9, 1.10, 2.4, 2.7)
Содержание
- 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 Общее описание алгоритма
Алгоритм Ланцоша был опубликовн физиком и математиком Корнелием Ланцощем в 1950 году[1]. Этот метод является частным случаем алгоритма Арнольда в случае, если исходная матрица [math]A[/math] - симметрична, и был представлен как итерационный метод вычисления собственных значений симметричной матрицы. Этот метод позволяет за [math]k[/math] итераций вычислять [math]k[/math] приближений собственных значений и собственных векторов исходной матрицы. Хотя алгоритм и был эффективным в вычислительном смысле, но он на некоторое время был предан забвению из-за численной неустойчивости. Только в 1970 Ojalvo и Newman модифицировали алгоритм для использования в арифметике с плавающей точкой[2]. Новый метод получил название алгоритма Ланцоша с полной переортогонализацией. Но эта статья про его исходную версию. На вход алгоритма подается вещественная симметричная матрица [math]A = A^{T}[/math],
- [math] A = \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1\ n-1} & a_{1\ n} \\ a_{12} & a_{22} & a_{23} & \cdots & a_{2\ n-1} & a_{2\ n} \\ a_{13} & a_{23} & a_{33} & \cdots & a_{3\ n-1} & a_{3\ n} \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ a_{1\ n-1} & \cdots & \cdots & a_{n-2\ n-1} & a_{n-1\ n-1} & a_{n-1\ n} \\ a_{1\ n} & \cdots & \cdots & a_{n-2\ n} & a_{n-1\ n} & a_{n\ n} \\ \end{pmatrix} [/math]
Поэтому достаточно хранить только чуть больше половины элементов исходной матрицы.
Сам алгоритм соединяет в себя метод Ланцоша построения крыловского подпространства с процедурой Релея-Ритца. На каждой итерации строится матрица [math]Q_k = [q_1, q_2, \dots, q_k][/math] размерности [math]n \times k[/math], состоящая из ортонормированных векторов Ланцоша. А в качестве приближенных собственных значений берутся числа Ритца, т.е. собственные значения симметричной трехдиагональной матрицы [math]T_k = Q^T_k A Q[/math] размерности [math]k \times k[/math].
- [math] T_k = \begin{pmatrix} \alpha_1 & \beta_1 \\ \beta_1 & \alpha_2 & \beta_2 \\ & \beta_2 & \ddots & \ddots \\ & & \ddots & \ddots & \beta_{k-1} \\ & & & \beta_{k-1} & \alpha_k \end{pmatrix} [/math]
Нахождение собственных значений и собственных векторов такой матрицы значительно проще чем их вычисление для исходной матрицы. Например, они могут быть вычислены с помощью метода «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы. В этом методе эта процедура занимает [math]O(k^3)[/math] операций, где константа оказывается на практике довольно мала.
1.2 Математическое описание алгоритма
Исходные данные: симметрическая матрица [math]A[/math], начальный вектор [math]b[/math].
Вычисляемые данные: собственные вектора матрицы [math]T_k[/math] являющиеся столбцами матрицы [math]Q_k V[/math], и матрица собственных значений [math]\Lambda[/math], где [math]V, \Lambda[/math] из спектрального разложения [math]T_k = V\Lambda V^T[/math]. До начала первой итерации нормируется вектор начального приближения [math]q_1 = \frac{b}{\|b\|_2}[/math]. Затем задаются константы [math]\beta_0 = 0,\; q_0 = 0 [/math].
После выполнения данных действий начинаются итерации алгоритма. Которых будет не больше чем [math]k[/math]. В ходе итераций постепенно формируется матрица [math]T_k[/math] из [math]\alpha_i[/math]-x и [math]\beta_i[/math]-х.
На [math]i[/math]-й итерации происходят следующие вычисления. Вначале вычисляется произведение матрицы на вектор [math]z = Aq_i[/math]. Затем считается скалярное произведение [math]\alpha_i = q^T_i z[/math] и значение сохраняется как элемент формируемой матрицы. После этого необходимо произвести вычисление линейной комбинации векторов [math]z = z - \alpha_i q_i - \beta_{i-1}q_{i-1}[/math]. Норму полученного вектора запишем как элемент матрицы [math]\beta_i = \|z\|_2[/math]. Если [math]\beta_i[/math] оказалась равной нулю, то больше никаких итераций в алгоритме не производят, сразу переходят к вычислению собственных значений и собственных векторов. Иначе, в конце итерации считают [math]q_{i+1} = \frac{z}{\beta_i}[/math] и переходят к следующей итерации.
После всех итераций необходимо вычислить собственные значения и собственные вектора матрицы [math]T_k[/math].[3]
1.3 Вычислительное ядро алгоритма
Вычислительным ядром на каждой итерации является вычисление произведения матрицы на вектор:
- [math]z = Aq_i[/math]
1.4 Макроструктура алгоритма
Для удобства понимания алгоритма можно выделить следующие макрооперации:
- умножение матрицы на вектор. Состоит из операций умножения вектора на число и сложения векторов.
- скалярное произведение векторов.
- линейная комбинация векторов. Состоит из умножений вектора на число и сложений векторов.
- получение нормы вектора. Состоит из скалярного произведения векторов, вычисления квадратного корня.
- деление вектора на число.
- вычисление собственных векторов и собственных значений трехдиагональной симметричной матрицы.
На каждой итерации будет вычисляться умножение матрицы [math]A[/math] размера [math]n \times n[/math] на вектор [math]q_i[/math]. Эта операция и будет самой ресурсозатратной на каждой итерации. Результаты этой макрооперации будут использованы затем в скалярном произведении векторов и в линейной комбинации из 3 векторов. После этого результат линейной комбинации векторов необходимо нормировать. Для этого можно разбить эту процедуру на вычисление нормы вектора и деление вектора на получившееся число. Кроме того, после завершения итераций необходимо вычислить собственные вектора и собственные значения полученной матрицы. Данная операция не конкретизируется алгоритмом, и может быть выполнена любым способом. Поэтому удобнее её представить в виде макрооперации. В разделе 1.7 можно наглядно увидеть как передаются данные между этими макрооперациями.
1.5 Схема реализации последовательного алгоритма
Последовательность исполнения метода следующая:
[math]1.\, \beta_0 = 0,\; q_0 = 0[/math] #Инициализируются константы
[math]2.\, \|b\|_2 = \sqrt{\sum\limits_{j=1}^{n} b_j^2}[/math] #Вычисляем норму вектора начального приближения.
[math]3.\, q_{1_{j}} = \frac{b_{j}}{\|b\|_2}, \; j = 1,\, \dots\,, n[/math] #Нормализуем вектор начального приближения.
Начинаются итерации цикла, которых не больше чем [math]k[/math]. На [math]i[/math]-й итерации производятся следующие вычисления:
[math]i.1\, z_j = \sum\limits_{m=1}^{n} a_{jm} q_{i_m}, \; j = 1,\,\dots\,, n[/math] #Считаем результат применения линейного оператора [math]A[/math] к вектору [math]q_i[/math].
[math]i.2\, \alpha_i = \sum\limits_{j=1}^{n}q_{i_j} z_j[/math] #Получаем результат скалярного произведения векторов [math]q_i[/math] и [math]z[/math].
[math]i.3\, z_j = z_j - \alpha_i q_{i_j} - \beta_{i-1}q_{i-1_j}, \, j = 1,\,\dots\,, n[/math] #Вычисляем линейную комбинацию векторов.
[math]i.4\, \beta_i = \|z\|_2 = \sqrt{\sum\limits_{j=1}^{n} z_j^2}[/math] #Считаем норму вектора [math]z[/math].
[math]i.5[/math] Проверка равенства [math]\beta_i == 0[/math] # Если норма оказалась равной нулю, то завершаем итерации и переходим к вычислению собственных векторов и собственных значений полученной матрицы. В обратном случае, продолжаем выполнения итераций.
[math]i.6\, q_{i+1_j} = \frac{z_j}{\beta_i}, \; j = 1,\, \dots \,, n[/math] #Нормируем вектор [math]z[/math].
[math]i.7\,[/math] Если выполнили [math]k[/math] итераций, то завершаем выполнение итераций, переходим к следующему шагу. Иначе начинаем последующую итерацию цикла.
[math]4.[/math] Вычисляем собственные значения и собственные вектора полученной матрицы [math]T_k[/math].
1.6 Последовательная сложность алгоритма
Вычисления [math]k[/math] собственных значений матрицы порядка [math]n[/math] и соответствующих им собственных векторов алгоритмом Ланцоша без переортогонализации состоит из инициализации, [math]k[/math] итераций и нахождения собственных значений и векторов трехдиагональной симметричной матрицы порядка [math]k[/math].
Инициализация состоит из следующих операций:
- вычисление нормы:
- [math]n[/math] операций умножения
- [math]n-1[/math] операций сложений
- извлечение квадратного корня
- [math]n[/math] операций делений (нормирование)
Каждая из [math]k[/math] итераций состоит из:
- умножение матрицы на вектор:
- [math]n^2[/math] умножений
- [math]n(n-1)[/math] сложений
- скалярное произведение двух векторов
- [math]n[/math] умножений
- [math]n-1[/math] сложение
- вычисление линейной комбинации векторов:
- [math]2n[/math] умножений
- [math]2n[/math] вычитаний
- норма вектора:
- [math]n[/math] операций умножения
- [math]n-1[/math] операций сложений
- извлечение квадратного корня
- нормирование вектора:
- [math]n[/math] операций деления
Итого (без учета решения задачи собственных значений трехдиагональной матрицы):
- [math] k(n^2 + 4n) + n [/math] умножений,
- [math] k(n^2+3n-2) + n-1[/math] сложений/вычитаний,
- [math] kn + n [/math] делений,
- [math] k + 1 [/math] вычислений квадратного корня.
Умножения и сложения (вычитания) составляют основную часть алгоритма.
Для нахождение собственных значений и векторов трехдиагональной симметричной матрицы эффективней всего использовать метод «разделяй и властвуй» вычисления собственных значений и векторов симметричной трехдиагональной матрицы, который требует [math] O(k^3) [/math] операций. При классификации по последовательной сложности, таким образом, метод Ланцоша без переортогонализации относится к алгоритмам с квадратичной сложностью.
1.7 Информационный граф
1.8 Ресурс параллелизма алгоритма
Алгоритм Ланцоша в параллельной форме состоит из инициализации, [math]k[/math] итераций и вычисления собственных значений матрицы [math]T_k[/math], рассмотрение которого выходит за рамки данной статьи.
Заметим, что вычисление суммы [math]n[/math] элементов имеет высоту [math]\log n[/math] и линейную ширину ярусов [math]\frac{n}{2}, \frac{n}{4}, ... , 1[/math].
Инициализации состоит из следующих операций:
- вычисление нормы:
- вычисление скалярного произведения вектора самого на себя:
- ярус умножений шириной [math]n[/math]
- [math]\log{n}[/math] ярусов сложений с наибольшей шириной [math]\frac{n}{2}[/math]
- ярус извлечения квадратного корня с единичным вычислением
- вычисление скалярного произведения вектора самого на себя:
- деление вектора на число (нормирование вектора):
- ярус [math]n[/math] операций деления
Итерационная часть алгоритма состоит из следующих операций:
- умножение матрицы на число:
- ярус умножений с шириной [math]n^2[/math]
- [math]\log n[/math] ярусов с наибольшей шириной [math]\frac{n^2}{2}[/math] (блок вычисления n сумм n элементов)
- вычисление скалярного произведения двух векторов:
- ярус умножений шириной [math]n[/math]
- [math]\log{n}[/math] ярусов сложений с наибольшей шириной [math]\frac{n}{2}[/math]
- вычисление линейной комбинации векторов:
- ярус умножений шириной [math]2n[/math]
- два яруса сложений шириной [math]n[/math]
- вычисление нормы:
- вычисление скалярного произведения вектора самого на себя:
- ярус умножений шириной [math]n[/math]
- [math]\log{n}[/math] ярусов сложений с максимальной шириной [math]\frac{n}{2}[/math]
- ярус извлечения квадратного корня с единичным вычислением
- вычисление скалярного произведения вектора самого на себя:
- проверка на выход из цикла
- деление вектора на число (нормирование вектора, может отсутствовать на последней итерации):
- ярус [math]n[/math] операций деления
Таким образом, в параллельном варианте основную долю времени будут занимать операции сложения при умножении матрицы на вектор.
При классификации по высоте ЯПФ Алгоритм Ланцоша относится к алгоритмам со сложностью [math]O(k \log n)[/math]. При классификации по ширине ЯПФ его сложность будет [math]O(n^2)[/math].
1.9 Входные и выходные данные алгоритма
Входные данные: симметричная матрица [math]A[/math] (элементы [math]a_{ij}, i \geq j[/math]), число [math]k[/math].
Объём входных данных: [math]\frac{n (n + 1)}{2} + 1[/math] (в силу симметричности достаточно хранить только диагональ и над/поддиагональные элементы, единица относится к параметру [math]k[/math]).
Выходные данные: по [math]k[/math] приближений собственных значений и собственных векторов.
Объём выходных данных: [math]k+kn[/math].
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является субквадратичным: [math]\frac{kn^2}{k \log{n}}[/math].
Вычислительная мощность алгоритма нахождения коэффициентов [math]\alpha_i,\beta_i[/math] без экономии памяти состовляет [math]\frac{k(2n^2+8n-1)+3n}{n^2+2k}[/math], при [math]k \ll n[/math], вычислительная мощность приблизительно равна [math]2k[/math].
В случае хранения половины исходной матрицы получаем соотношношение [math]\frac{k(2n^2+8n-1)+3n}{n(n+1)/2+2k}[/math], что приблизительно равно [math]4k[/math].
Алгоритм Ланцоша без переортогонализации не является детерминированным (возможно выполнение меньшего числа итераций алгоритма), в случае если все собственные значения уже вычислены.
Важное свойство метода Ланцоша состоит в том, что первыми в матрице [math]T_{j}[/math] появляются собственные значения с максимальной величиной по модулю. Таким образом, метод особенно хорошо подходит для вычисления собственных значений матрицы [math]A[/math], находящихся на краях её спектра.
2 Программная реализация
2.1 Особенности реализации последовательного алгоритма
2.2 Локальность данных и вычислений
2.3 Возможные способы и особенности параллельной реализации алгоритма
2.4 Масштабируемость алгоритма и его реализации
Проведём исследование масштабируемости параллельной реализации алгоритма Ланцоша без переортогонализации согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета.
2.4.1 Реализация на языке C
Исследованная параллельная реализация на языке C
2.4.2 Масштабируемость реализации алгоритма
Сборка осуществлялась со следующими параметрами:
- gcc-5.2.0
- openmpi-1.8.4
- аргументы компилятора: -std=c11 -Ofast
Набор и границы значений изменяемых параметров реализации алгоритма:
- число процессоров [32 : 256] с шагом 16;
- размер матрицы [5000 : 200000] с шагом 5000.
В результате проведённых экспериментов был получен следующий диапазон эффективности реализации алгоритма:
- минимальная эффективность реализации 1.92%;
- максимальная эффективность реализации 59.47%.
График полученного распределения эффективности:
На следующих рисунках приведены графики производительности и эффективности данной реализации в зависимости от изменяемых параметров запуска.
Средняя эффективность для матриц с размером больше 25000 составила 52.7%.
2.5 Динамические характеристики и эффективность реализации алгоритма
2.6 Выводы для классов архитектур
2.7 Существующие реализации
Реализация в проекте IETL[4]
Реализация в проекте ARPACK [5]
3 Ссылки
- ↑ Lanczos, C. "An iteration method for the solution of the eigenvalue problem of linear differential and integral operators", J. Res. Nat’l Bur. Std. 45, 255-282 (1950).
- ↑ Ojalvo, I.U. and Newman, M., "Vibration modes of large structures by an automatic matrix-reduction method", AIAA J., 8 (7), 1234–1239 (1970).
- ↑ Деммель Д. Вычислительная линейная алгебра
- ↑ http://www.comp-phys.org/software/ietl/lanczos.html
- ↑ http://www.caam.rice.edu/software/ARPACK/