Текущая версия на 11:09, 19 декабря 2016

Основные авторы описания: А.А.Ульянов (разделы 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.9, 2.7), К.Петров (разделы 1.7, 1.8, 1.10, 2.4).

Обобщенный метод минимальных невязок
Последовательный алгоритм
Последовательная сложность	[math] O(nm^2) [/math]
Объём входных данных	[math] (n + 1)^2 [/math]
Объём выходных данных	[math] n [/math]

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

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

Обобщённый метод минимальных невязок (англ. Generalized minimal residual method, GMRES) - итерационный метод численного решения системы линейных алгебраических уравнений с невырожденной матрицей. Метод основан на минимизации квадратичного функционала невязки на подпространствах Крылова. Разработан Юсефом Саадом и Мартином Шульцем^[1] в 1986 году как обобщение метода MINRES^[2] на случай систем с несимметричными матрицами.

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

Исходные данные:

система линейных алгебраических уравнений вида [math] Ax = b [/math], где [math] A [/math] — невырожденная матрица размера [math]n[/math]-на-[math] n [/math];
[math] x_0 [/math] - начальное приближение.

Вычисляемые данные:

[math] x_m [/math] - приближённое решение исходной системы.

Метод GMRES приближает точное решение исходной системы [math] Ax = b [/math] вектором [math] x_m [/math], который минимизирует Евклидову норму невязки [math] r_m = Ax_m - b[/math] на подпространстве Крылова:

[math] K_m = K_m(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{m-1}b \} [/math], где [math] \operatorname{span} \, \{ a_1, a_2, \ldots, a_m \} [/math] - линейная оболочка системы векторов [math] \{ a_1, a_2, \ldots, a_m \} [/math].

На каждой итерации метода решение уточняется поправкой, представленной в виде разложения по ортонормированному базису пространства [math] K_m [/math]:

[math] x_m = x_0 + z_m [/math], где [math] x_0 [/math] - некоторое начальное приближение, [math] z_m \in K_m [/math] - поправка решения.

1.2.1 Ортогонализация Арнольди

Для построения ортонормированного базиса [math] v_1, v_2, \ldots, v_m [/math] пространства Крылова [math] K_m [/math] метод использует процесс Арнольди.

Идея процесса Арнольди состоит в том, что для нахождения вектора [math] v_{j+1} [/math] при известных векторах [math] v_1, v_2, \ldots, v_j [/math] необходимо вычислить вектор [math] Av_j [/math] и ортонормировать этот вектор относительно [math] v_1, v_2, \ldots, v_j [/math]:

[math] v_{j+1} = \frac{\hat{v}_{j+1}}{\|\hat{v}_{j+1}\|_2} [/math], где [math] \hat{v}_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{i} [/math].

Коэффициенты [math] h_{ij} = (Av_j, v_i) [/math] определяются из условия ортогональности:

[math] A(v_j, v_k) - \sum_{i=1}^j h_{ij}(v_{i}, v_k) = 0 [/math], где [math] \quad k=1,\ldots,j [/math].

Результат процедуры Арнольди совпадает с результатом процесса Грамма-Шмидта, применённого к векторам [math] \{b, Ab, A^2b, \ldots, A^{m-1}b\} [/math].

Исходные данные:

[math] v_1 [/math], такой что [math] \|v_1\|_2 = 1 [/math];
матрица [math] A [/math] размером [math] n [/math]-на-[math] n [/math].

Вычисляемые данные:

[math] v_1, v_2, \ldots, v_m [/math] - ортонормированный базис пространства Крылова [math] K_m [/math].

Алгоритм процесса:

Выполнять для [math] \quad j=1,\ldots,m [/math]:

[math] h_{ij} = (Av_j, v_i), \quad i=1,\ldots,j [/math];
[math] \hat{v}_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{i} [/math];
[math] h_{j+1j} = \|\hat{v}_{j+1}\|_2 [/math];
[math] v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1j}} [/math].

При введении матричных обозначений можно записать:

Поправка решения:

[math] z_m = V_my_m [/math], где [math] y_m \in \mathbb{R}^m [/math] - вектор коэффициентов, [math] V_m = [v_1|v_2|...|v_m] [/math].

Матричное представление процесса Арнольди:

[math] AV_m = V_{m+1}H_{m+1m} = V_mH_m + h_{m+1m}v_{m+1} e_m^T [/math],

где

[math]H_{m+1m} = \begin{bmatrix}h_{11} & h_{12} & h_{13} &... & h_{1m-1} & h_{1m} \\ h_{21} & h_{22} & h_{23} & ... & h_{2m-1} & h_{2m} \\ 0 & h_{32} & h_{33} &...& h_{3m-1} & h_{3m} \\ ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & ... & h_{m m-1} & h_{mm} \\0 & 0 & 0 & ... & 0 & h_{m+1 m}\end{bmatrix}[/math]

- верхняя матрица Хессенберга размерности [math] m+1 [/math] на [math] m [/math].

1.2.2 Минимизация функционала невязки

Невязка исходной системы имеет вид:

[math] r_m = b - A(x_0 + z_m) = r_0 - AV_my_m = \beta v_1 - V_{m+1}H_{m+1m}y_m[/math], где [math] \beta = \|r_0\|_2[/math].

C учетом того, что матрица [math] V_{k+1} [/math] составлена из ортонормированных векторов:

[math] r_m = \beta V_{m+1}e_1 - H_{m+1m}y_m[/math], где [math] e_1 = (1,0,0,\ldots,0)^T [/math].

Норма невязки [math] r_m [/math]:

[math] \|r_m\|_2 = \|\beta e_1 - H_{m+1m}y_m\|_2 [/math].

Норма [math] \|r_m\|_2 [/math] достигает минимума при [math] \beta e_1 = H_{m+1m}y_m [/math].

Для нахождения решения такой системы, матрица [math] H_{m+1m} [/math] приводится к верхнему треугольному виду с помощью [math] m [/math] последовательных вращений Гивенса.

1.2.3 Общая схема метода

В общем виде m-aя итерация алгоритма GMRES может быть записан следующим образом:

найти ортонормированный базис [math] V_m [/math] подпространства [math] K_m [/math] с помощью ортогонализации Арнольди;
найти [math] y_m [/math], минимизирующий [math] \|r_m\|_2 [/math];
вычислить [math] x_m = x_0 + V_my_m [/math], где [math] x_0 [/math] - начальное приближение;
вычислить [math] r_m [/math] и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода GMRES состоит из двух частей:

Вычисление ортонормированного базиса [math] K_m [/math];
Формирование приближенного решения [math] x_m [/math].

Для вычисления ортонормированного базиса [math] K_m [/math] метод использует процесс Арнольди:

[math] h_{j+1j}v_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{j} [/math], где [math] h_{ij} = (Av_j, v_i)[/math], [math] \quad j=1,\ldots,m [/math].

Для нахождения приближённого решения [math] x_m [/math] метод использует формулу:

[math] x_m = x_0 + V_my_m [/math].

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

В алгоритме можно выделить следующие макрооперации:

Умножение матрицы на вектор;
Вычисление скалярного произведения векторов;
Вычисление Евклидовой нормы вектора;
Умножение вектора на скаляр.

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

На практике используется перезапускаемая версия метода - GMRES(m):

1 Подготовка перед итерационным процессом:

1.1 Выбрать начальное приближение [math] x_0 [/math];

1.2 Посчитать невязку [math] r_0 = b - Ax_0 [/math];

1.3 Вычислить [math] v_1 = \frac{r_0}{\|r_0\|_2} [/math].

2 Построение ортонормированного базиса [math] K_m [/math]:

Для всех [math] j [/math] от 1 до m по нарастанию выполнять:

2.1 [math] h_{ij} := (Av_j, v_i), \quad i=1,\ldots,j [/math];

2.2 [math] \hat{v}_{j+1} := Av_j - \sum_{i=1}^j h_{ij}v_{i} [/math];

2.3 [math] h_{j+1j} = \|\hat{v}_{j+1}\|_2 [/math];

2.4 [math] v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1j}} [/math].

3 Вычисление [math] x_m [/math]:

3.1 [math] x_m = x_0 + V_my_m [/math], где [math]y_m[/math] минимизирует [math]\|r_0 - AV_my_m\|_2[/math];

3.2 Вычислить [math] r_m [/math];

3.3 Если требуемая точность достигнута, остановиться.

4 Перезапуск:

4.1 [math]x_0 = x_m[/math];

4.2 [math]v_1 = \frac{r_m}{\|r_m\|_2}[/math];

4.3 Перейти к шагу 2.

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

Если пренебречь сложностью вычисления [math] y_m [/math], то общую сложность алгоритма GMRES можно разделить на две части:

1. Сложность вычисления ортонормированного базиса пространства [math] K_m [/math]:

a) Для вычисления [math] j [/math]-го вектора базиса [math] K_m, j \lt m [/math] требуется:

[math] NZ + (2j + 1)n [/math] мультипликативных операций, где [math] NZ [/math] - количество ненулевых элементов матрицы [math] A [/math].

b) Вычисление последнего вектора базиса требует:

[math] n(m + 1) [/math] мультипликативных операций.

c) Общая мультипликативная сложность вычисления ортонормированного базиса [math] K_m [/math]:

[math] O(nm^2) [/math].

2. Сложность вычисления приближённого решения [math] x_m [/math]:

a) Вычисление этой формулы требует [math] nm [/math] мультипликативных операций.

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

Рис. 1 input - входные данные алгоритма (квадратная матрица, вектор правой части, начальное приближение, точность решения); [math] A_i [/math] - шаги алгоритма Арнольди, [math] x_i [/math] - нахождение невязки на [math] i [/math]-ой итерации; if - вычисление [math] r_m [/math] и в случае достижения нужной точности - остановка.

Шаги алгоритма Арнольди выполняются строго последовательно, но вычисление [math] i [/math]-ой невязки может происходить параллельно с [math] (i + 1) [/math]-ым шагом алгоритма Арнольди.

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

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

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

При выполнении параллельного алгоритма решения систем линейных алгебраических уравнений необходимо распределить входные данные наиболее рациональным способом. Матрица коэффициентов разбивается на непрерывные горизонтальные полосы. Распределение вектора правой части и начального приближения выполнено по принципу дублирования, т.е. все элементы векторов скопированы на все процессоры. Для рассматриваемой задачи решения СЛАУ такой подход распределения допустим, т.к. векторы содержат столько же элементов, сколько и одна строка или один столбец матрицы коэффициентов.

Подзадачи можно разделить на 2 группы: операции между матрицей и вектором (например, умножение матрицы на вектор) и операции между двумя векторами (например, скалярное умножение). Подзадачи умножения константы на вектор и определение нормы вектора можно относить ко второй группе, так как их трудоемкость одного порядка.

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

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

Входные данные:

[math] A [/math] – квадратная матрица размера [math] n [/math]-на-[math] n [/math];
[math] b [/math] – вектор правой части размера [math] n [/math];
[math] x_0 [/math] – начальное приближение размера [math] n [/math];
[math] \epsilon [/math] – погрешность решения.

Для перезапускаемой версии алгоритма дополнительно:

[math] m [/math] – число итераций.

Объём входных данных: [math] n^2 + 2n + 1 [/math].

Выходные данные:

[math] x_m [/math] - вектор приближённого решения размера [math] n [/math].

Объём выходных данных: [math] n [/math].

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. GMRES использует перезазапуски для контроля требований к хранению данных.

Если используется неперезапускаемый метод, GMRES будет сходиться за не более [math] n [/math] шагов. Но это не имеет практического применения, когда [math] n [/math] велико. Кроме того, требования к памяти и вычислениям при отсутствии перезапусков являются очень высокими. Важнейшим элементом для успешного применения [math] GMRES(m)[/math] являются перезапуски; то есть выбор [math] m [/math]. Существуют примеры, для которых метод застаивается и сходимость есть только на [math]m[/math]-м шаге. Для таких систем, любой выбор [math] m [/math] меньше, чем [math] n [/math] - не сходится.

Основным недостатком GMRES является то, что объемы работ и ресурсы памяти, требуемые для итерации, линейно возрастают с увеличением числа итераций. Если не удалось получить чрезвычайно быструю сходимость - затраты будут очень высокими. Основной способ преодолеть это ограничение - это перезапуск итерации. После выбранного числа итераций [math] m [/math] накопленные данные будут очищены и промежуточные результаты используются в качестве входных данных для последующих итераций. Эта процедура повторяется до тех пор, пока сходимость не будет достигнута. Трудность заключается в выборе подходящего значения [math] m [/math]. Если [math] m [/math] слишком мало, [math]GMRES(m)[/math] может быть медленно сходиться, или не сойтись совсем. Значение [math] m [/math], большее, чем это необходимо, включает в себя чрезмерную работу и использует больше памяти. К сожалению, нет определенных правил, регулирующих выбор [math] m [/math]; выбирая , когда делать перезапуск - это вопрос опыта.

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

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

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

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

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

Масштабируемость рассматривалась в работах ^[3] и ^[4]. Для исследования масштабируемости была использована параллельная реализация с использованием объекта KSP(решатель СЛАУ) бибилиотеки PETSc(содержащей структуры, типы данных и подпрограммы для параллельной реализации методов на подпространствах Крылова), найденная на GitHub ^[5].

Вычисления проводились на вычислительном кластере, состоящем из 16 узлов, со следующей конфигурацией: Altus 1804i Server - 4P Interlagos Node, Quad AMD Opteron 6276, 16C, 2.3 GHz, 128GB, DDR3-1333 ECC, 80GB SSD, MLC, 2.5" HCA, Mellanox ConnectX 2, 1-port QSFP, QDR, memfree, CentOS, Version 5.

На рис. 2 представлен график сильной масштабируемости.

Рис. 2 Зависимость производительности от числа процессоров (запуск производился на 1, 8, 27, 64, 216, 512, 1000 процессорах)

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:

3 Литература

↑ Y.Saad, M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Scientific and Stat. Comp. 7: 856-869 (1986)
↑ C. C. Paige and M. A. Saunders. Solution of sparse indefinite systems of linear equations, SIAM J. Numerical Analysis 12, 617-629 (1975)
↑ https://arxiv.org/pdf/1606.08740.pdf
↑ https://smartech.gatech.edu/bitstream/handle/1853/55635/PRATAPA-DISSERTATION-2016.pdf
↑ https://github.com/erdc-cm/petsc-dev/tree/master/src/ksp/ksp/impls/gmres

[1] Y.Saad, M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Scientific and Stat. Comp. 7: 856-869 (1986)

[2] C. C. Paige and M. A. Saunders. Solution of sparse indefinite systems of linear equations, SIAM J. Numerical Analysis 12, 617-629 (1975)

[3] ttps://arxiv.org/pdf/1606.08740.pdf

[4] ttps://smartech.gatech.edu/bitstream/handle/1853/55635/PRATAPA-DISSERTATION-2016.pdf

[5] ttps://github.com/erdc-cm/petsc-dev/tree/master/src/ksp/ksp/impls/gmres

[1]

[2]

[3]

[4]

[5]

@@ Строка 1: / Строка 1: @@
-{{Assignment}}
 Основные авторы описания: [[#Участник:Alexander34396|А.А.Ульянов]] (разделы [[#Общее описание алгоритма|1.1]], [[#Математическое описание алгоритма|1.2]], [[#Вычислительное ядро алгоритма|1.3]], [[#Макроструктура алгоритма|1.4]], [[#Схема реализации последовательного алгоритма|1.5]], [[#Последовательная сложность алгоритма|1.6]], [[#Входные и выходные данные алгоритма|1.9]], [[#Существующие реализации алгоритма|2.7]]), [[#Участник:Kpetrov|К.Петров]] (разделы [[#Информационный граф|1.7]], [[#Ресурс параллелизма алгоритма|1.8]], [[#Свойства алгоритма|1.10]], [[#Масштабируемость алгоритма и его реализации|2.4]]).
@@ Строка 9: / Строка 9: @@
 }}
-== Свойства и структура алгоритма ==
+= Свойства и структура алгоритмов =
-=== Общее описание алгоритма ===
+== Общее описание алгоритма ==
 '''Обобщённый метод минимальных невязок''' ({{lang-en|Generalized minimal residual method, GMRES}}) - итерационный метод численного решения системы линейных алгебраических уравнений с невырожденной матрицей. Метод основан на минимизации квадратичного функционала невязки на подпространствах Крылова. Разработан Юсефом Саадом и Мартином Шульцем<ref> Y.Saad, M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Scientific and Stat. Comp. 7: 856-869 (1986) </ref> в 1986 году как обобщение метода MINRES<ref>C. C. Paige and M. A. Saunders. Solution of sparse indefinite systems of linear equations, SIAM J. Numerical Analysis 12, 617-629 (1975)</ref> на случай систем с несимметричными матрицами.
-=== Математическое описание алгоритма ===
+== Математическое описание алгоритма ==
 '''Исходные данные:'''
@@ Строка 23: / Строка 23: @@
 * <math> x_m </math> - приближённое решение исходной системы.
-Метод GMRES приближает точное решение исходной системы <math> Ax = b </math> вектором <math> x_m </math>, который минимизирует <math> l_2 </math>-норму невязки <math> r_m = Ax_m - b</math> на подпространстве Крылова:
+Метод GMRES приближает точное решение исходной системы <math> Ax = b </math> вектором <math> x_m </math>, который минимизирует Евклидову норму невязки <math> r_m = Ax_m - b</math> на подпространстве Крылова:
-: <math> K_m = K_m(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{m-1}b \} </math>.
+: <math> K_m = K_m(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{m-1}b \} </math>, где <math> \operatorname{span} \, \{ a_1, a_2, \ldots, a_m \} </math> - линейная оболочка системы векторов <math> \{ a_1, a_2, \ldots, a_m \} </math>.
-На каждой итерации метода решение уточняется поправкой, представленной в виде разложения по <math> l_2 </math>-ортонормированному базису пространства <math> K_m </math>:
+На каждой итерации метода решение уточняется поправкой, представленной в виде разложения по ортонормированному базису пространства <math> K_m </math>:
 : <math> x_m = x_0 + z_m </math>, где <math> x_0 </math> - некоторое начальное приближение, <math> z_m \in K_m </math> - поправка решения.
-==== Ортогонализация Арнольди ====
+=== Ортогонализация Арнольди ===
+Для построения ортонормированного базиса <math> v_1, v_2, \ldots, v_m </math> пространства Крылова <math> K_m </math> метод использует процесс Арнольди.
+Идея процесса Арнольди состоит в том, что для нахождения вектора <math> v_{j+1} </math> при известных векторах <math> v_1, v_2, \ldots, v_j </math> необходимо вычислить вектор <math> Av_j </math> и ортонормировать этот вектор относительно <math> v_1, v_2, \ldots, v_j </math>:
+: <math> v_{j+1} = \frac{\hat{v}_{j+1}}{\|\hat{v}_{j+1}\|_2} </math>, где <math> \hat{v}_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{i} </math>.
-Для построения ортонормированного базиса <math> K_m </math> метод использует процесс Арнольди.
+Коэффициенты <math> h_{ij} = (Av_j, v_i) </math> определяются из условия ортогональности:
+: <math> A(v_j, v_k) - \sum_{i=1}^j h_{ij}(v_{i}, v_k)  = 0 </math>, где <math> \quad k=1,\ldots,j </math>.
+Результат процедуры Арнольди совпадает с результатом процесса Грамма-Шмидта, применённого к векторам <math> \{b, Ab, A^2b, \ldots, A^{m-1}b\} </math>.
 '''Исходные данные:'''
@@ Строка 37: / Строка 45: @@
 * матрица <math> A </math> размером <math> n </math>-на-<math> n </math>.
-'''Формула процесса:'''
+'''Вычисляемые данные:'''
-: <math> h_{j+1j}v_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{j} </math>, где <math> h_{ij} = (Av_j, v_i)</math>, <math> \quad j=1,\ldots,m </math>.
+* <math> v_1, v_2, \ldots, v_m </math> - ортонормированный базис пространства Крылова <math> K_m </math>.
 '''Алгоритм процесса:'''
@@ Строка 49: / Строка 57: @@
 '''При введении матричных обозначений можно записать:'''
-# <math> z_m = V_my_m </math>, где <math> y_m \in \mathbb{R}^m </math> - вектор коэффициентов;
-# <math> AV_m = V_mH_m + h_{m+1m}v_{m+1} e_m^T = V_{m+1}\overline{H}_m </math>, где <math> V_m = [v_1|v_2|...|v_m] </math>, а <math> \overline{H}_m </math> - верхняя матрица Хессенберга размерности <math> (j+1) </math> на <math> j </math>.
-==== Минимизация функционала невязки ====
+Поправка решения:
+: <math> z_m = V_my_m </math>, где <math> y_m \in \mathbb{R}^m </math> - вектор коэффициентов, <math> V_m = [v_1|v_2|...|v_m] </math>.
+Матричное представление процесса Арнольди:
+: <math> AV_m = V_{m+1}H_{m+1m} = V_mH_m + h_{m+1m}v_{m+1} e_m^T </math>,
+где
+: <math>H_{m+1m} = \begin{bmatrix}h_{11} & h_{12} & h_{13} &... & h_{1m-1} & h_{1m} \\ h_{21} & h_{22} & h_{23} & ... & h_{2m-1} & h_{2m} \\ 0 & h_{32} & h_{33} &...& h_{3m-1} & h_{3m} \\ ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & ... & h_{m m-1} & h_{mm} \\0 & 0 & 0 & ... & 0 & h_{m+1 m}\end{bmatrix}</math>
+- верхняя матрица Хессенберга размерности <math> m+1 </math> на <math> m </math>.
+=== Минимизация функционала невязки ===
+Невязка исходной системы имеет вид:
+: <math> r_m = b - A(x_0 + z_m) = r_0 - AV_my_m = \beta v_1 - V_{m+1}H_{m+1m}y_m</math>, где <math> \beta = \|r_0\|_2</math>.
+C учетом того, что матрица <math> V_{k+1} </math> составлена из ортонормированных векторов:
+: <math> r_m = \beta V_{m+1}e_1 - H_{m+1m}y_m</math>, где <math> e_1 = (1,0,0,\ldots,0)^T </math>.
-Для решения исходной системы GMRES вычисляет приближённое решение <math> x_m </math>, которое минимизизует функционал невязки:
+Норма невязки <math> r_m </math>:
-: <math> J(y) = \|b- A x_m\|_2 = \|b - A(x_0 + V_my)\|_2 = \|\beta e_1 - \overline{H}_my_m \|_2 </math>, где <math> \beta = \|r_0\|_2 </math>.
+: <math> \|r_m\|_2 = \|\beta e_1 - H_{m+1m}y_m\|_2 </math>.
-Вектор <math> y_m </math> может найден как решение линейной задачи наименьших квадратов размера <math> m+1 </math>-на-<math> m </math>:
+Норма <math> \|r_m\|_2 </math> достигает минимума при <math> \beta e_1 = H_{m+1m}y_m </math>.
-: <math> y_m = \arg\min_{y}\| \beta e_1 - \overline{H}_my\|_2 </math>.
-Для решения задачи матрица <math> \overline{H}_m </math> приводится к верхнему треугольному виду с помощью <math> m </math> последовательных вращений Гивенса.
+Для нахождения решения такой системы, матрица <math> H_{m+1m} </math> приводится к верхнему треугольному виду с помощью <math> m </math> последовательных [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|вращений Гивенса]].
-==== Общая схема метода ====
+=== Общая схема метода ===
 В общем виде m-aя итерация алгоритма GMRES может быть записан следующим образом:
@@ Строка 71: / Строка 91: @@
 # вычислить <math> r_m </math> и остановиться,если требуемая точность была достигнута, иначе повторить.
-=== Вычислительное ядро алгоритма ===
+== Вычислительное ядро алгоритма ==
 Вычислительное ядро последовательной версии метода GMRES состоит из двух частей:
@@ Строка 83: / Строка 103: @@
 : <math> x_m = x_0 + V_my_m </math>.
-=== Макроструктура алгоритма ===
+== Макроструктура алгоритма ==
 В алгоритме можно выделить следующие макрооперации:
@@ Строка 89: / Строка 109: @@
 * Умножение матрицы на вектор;
 * Вычисление скалярного произведения векторов;
-* Вычисление <math> l_2 </math>-нормы вектора;
+* Вычисление Евклидовой нормы вектора;
-* Умножение вектора на скаляр;
+* Умножение вектора на скаляр.
-* Деление вектора на скаляр.
-=== Схема реализации последовательного алгоритма ===
+== Схема реализации последовательного алгоритма ==
 На практике используется перезапускаемая версия метода - GMRES(m):
@@ Строка 120: / Строка 139: @@
 :4.3 Перейти к шагу 2.
-=== Последовательная сложность алгоритма ===
+== Последовательная сложность алгоритма ==
 Если пренебречь сложностью вычисления <math> y_m </math>, то общую сложность алгоритма GMRES можно разделить на две части:
@@ Строка 136: / Строка 155: @@
-=== Информационный граф ===
+== Информационный граф ==
-----
-[[File: Drawing.png|300px|caption]]
-----
-На информационном графе:
-* '''input''' - входные данные алгоритма (квадратная матрица, вектор правой части, начальное приближение, точность решения);
-* <math> A_i </math> - шаги алгоритма Арнольди, <math> x_i </math> - нахождение невязки на <math> i </math>-ой итерации;
-* '''if''' - вычисление <math> r_m </math> и в случае достижения нужной точности - остановка.
-При этом, шаги алгоритма Арнольди выполняются строго последовательно, но вычисление <math> i </math>-ой невязки может происходить параллельно
+[[File: Drawing.png|thumb|center|350px|'''Рис. 1'''
+'''input''' - входные данные алгоритма (квадратная матрица, вектор правой части, начальное приближение, точность решения);
+ '''<math> A_i </math>''' - шаги алгоритма Арнольди, <math> x_i </math> - нахождение невязки на <math> i </math>-ой итерации;
+'''if''' - вычисление <math> r_m </math> и в случае достижения нужной точности - остановка.]]
+Шаги алгоритма Арнольди выполняются строго последовательно, но вычисление <math> i </math>-ой невязки может происходить параллельно
 с <math> (i + 1) </math>-ым шагом алгоритма Арнольди.
-=== Ресурс параллелизма алгоритма ===
+== Ресурс параллелизма алгоритма ==
 В ходе анализа перезапускаемого метода обобщенных минимальных невязок был выявлен ряд подзадач, выполняющихся неоднократно: умножение матрицы на вектор, вычисление скалярного произведения векторов, нахождение нормы вектора, нормирование вектора, нахождение взвешенной суммы векторов. Наиболее трудоемкой среди перечисленных операций является умножение матрицы на вектор.
-Выполнение перечисленных операций занимает большую часть времени работы алгоритма GMRES. Поэтому при увеличении эффективности использования вычислительных ресурсов и ресурсов памяти вычислительной машины при выполнении этих операций повышается эффективность метода в целом. Достигается это применением при выполнении перечисленных операций аппарата параллельных вычислений.
+Выполнение перечисленных операций занимает большую часть времени работы алгоритма GMRES. Поэтому здесь применяется аппарат параллельных вычислений для более эффективного использования ресурсов памяти вычислительной машины и для повышения эффективности метода в целом.
 При выполнении параллельного алгоритма решения систем линейных алгебраических уравнений необходимо распределить входные данные наиболее рациональным способом. Матрица коэффициентов разбивается на непрерывные горизонтальные полосы. Распределение вектора правой части и начального приближения выполнено по принципу дублирования, т.е. все элементы векторов скопированы на все процессоры. Для рассматриваемой задачи решения СЛАУ такой подход распределения допустим, т.к. векторы содержат столько же элементов, сколько и одна строка или один столбец матрицы коэффициентов.
@@ Строка 158: / Строка 174: @@
 Подзадачи можно разделить на 2 группы: операции между матрицей и вектором (например, умножение матрицы на вектор) и операции между двумя векторами (например, скалярное умножение). Подзадачи умножения константы на вектор и определение нормы вектора можно относить ко второй группе, так как их трудоемкость одного порядка.
-Результатом выполнения рассматриваемых операций является, вектор, части которого расположены на участвовавших в вычислениях процессорах. Для объединения результатов расчета и получения полного вектора на каждом процессоре вычислительной системы необходимо выполнить операцию сбора данных.
+Результатом выполнения рассматриваемых операций является вектор, части которого расположены на участвовавших в вычислениях процессорах. Для объединения результатов расчета и получения полного вектора на каждом процессоре вычислительной системы необходимо выполнить операцию сбора данных.
-=== Входные и выходные данные алгоритма ===
+== Входные и выходные данные алгоритма ==
 '''Входные данные:'''
@@ Строка 178: / Строка 194: @@
 '''Объём выходных данных''': <math> n </math>.
-=== Свойства алгоритма ===
+== Свойства алгоритма ==
+Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. GMRES использует перезазапуски для контроля требований к хранению данных.
+Если используется неперезапускаемый метод, GMRES будет сходиться за не более <math> n </math> шагов. Но это не имеет практического применения, когда <math> n </math> велико. Кроме того, требования к памяти и вычислениям при отсутствии перезапусков являются очень высокими. Важнейшим элементом для успешного применения <math> GMRES(m)</math> являются перезапуски; то есть выбор <math> m </math>. Существуют примеры, для которых метод застаивается и сходимость есть только на <math>m</math>-м шаге. Для таких систем, любой выбор <math> m </math> меньше, чем <math> n </math> - не сходится.
+Основным недостатком GMRES является то, что объемы работ и ресурсы памяти, требуемые для итерации, линейно возрастают с увеличением числа итераций. Если не удалось получить чрезвычайно быструю сходимость - затраты будут очень высокими. Основной способ преодолеть это ограничение - это перезапуск итерации. После выбранного числа итераций <math> m </math> накопленные данные будут очищены и промежуточные результаты используются в качестве входных данных для последующих итераций. Эта процедура повторяется до тех пор, пока сходимость не будет достигнута. Трудность заключается в выборе подходящего значения <math> m </math>. Если <math> m </math> слишком мало, <math>GMRES(m)</math>  может быть медленно сходиться, или не сойтись совсем. Значение <math> m </math>, большее, чем это необходимо, включает в себя чрезмерную работу и использует больше памяти. К сожалению, нет определенных правил, регулирующих выбор <math> m </math>; выбирая , когда делать перезапуск - это вопрос опыта.
-Вычислительная схема требует работы со следующими объектами: разреженной матрицей исходной СЛАУ — <math>n + 2k </math> вещественных и <math> n + k + 1 </math> целых (причем обычно двойной длины!) чисел; векторами правой части, искомого решения и невязки — <math>3n</math> вещественных чисел; векторами базиса подпространства Крылова — <math> mn </math> вещественных чисел; плотной матрицей малой переопределенной СЛАУ — <math> m(m + 1) </math> вещественных чисел; коэффициентами вращений Гивенса — <math> 2m </math> вещественных чисел; вектором решения малой СЛАУ (как с коэффициентами линейного комбинирования векторов базиса подпространства Крылова) — <math> m </math> вещественных чисел.
+= Программная реализация алгоритма =
-Полагая объем памяти для счетчиков циклов и временных переменных скалярного типа пренебрежимо малым, можно получить, что всего требуется хранить <math> n + k + 1 </math> целых чисел двойной длины и <math>4n + 2k + mn + (m + 1)2 + 3m + 1</math> — вещественных. При использовании обычной точности это требует <math>20n+12k+4mn+4(m+1)2+12m+8</math> байт, двойной — <math>36n + 20k + 8mn + 8(m + 1)2 + 24m + 16</math> байт. С учетом того, что <math>1 < m \ll n </math> и <math> k \ll n </math>, можно говорить, что требования метода GMRES к объему доступной памяти составляют <math> O(mn) </math>.
+== Особенности реализации последовательного алгоритма ==
-Основные вычислительные затраты определяются операциями нахождения скалярных произведений и линейного комбинирования векторов. Эти действия требуют <math> 2n^3 + \frac {3}{2}n^2- \frac{1}{2}n </math> арифметических операций. Доля соответствующих вычислительных затрат в схеме метода, таким образом, равна <math> \frac {2n^3 + \frac {3}{2}n^2- \frac{1}{2}n}{2n^3 + 7n^2 + 3n} </math>, т. е. при больших размерностях СЛАУ близка к 1.
+== Локальность данных и вычислений ==
-При этом строки матрицы СЛАУ, начиная с первой, последовательно исключаются из процесса решения. В такой ситуации наиболее эффективной является параллелизация, основанная на столбцовой декомпозиции матрицы СЛАУ при реализации вычислений на топологии межпроцессорных соединений “каждый с каждым”. При этом работу с вектором правой части целесообразно дублировать на каждом из процессоров. При использовании такой схемы на <math>p</math> процессорах число арифметических действий, выполняемых одним процессором, сокращается до величины <math> N^A_p(n) = \frac {n(3n^2+2n+np+p+1)}{p} </math>.
+== Возможные способы и особенности параллельной реализации алгоритма ==
-В вычислительной схеме GMRES присутствуют два интенсивных в вычислительном плане этапа: умножение матрицы СЛАУ на вектор (для нахождения невязок) и нахождение скалярных произведений векторов (при построении ортонормального базиса подпространства Крылова). Кроме того, именно на хранение векторов базиса <math> \left \{ v_j \right \} </math> в GMRES приходится бОльшая часть затрат памяти ( <math> mn </math> вещественных чисел, что часто оказывается больше объема памяти под разреженную матрицу СЛАУ). При этом в GMRES-схеме присутствует достаточно самостоятельная подзадача — решение переопределенной СЛАУ малой размерности в смысле наименьших квадратов, которая должна рассматриваться отдельно.
+== Масштабируемость алгоритма и его реализации ==
-== Программная реализация алгоритма ==
-=== Особенности реализации последовательного алгоритма ===
+Масштабируемость рассматривалась в работах <ref>https://arxiv.org/pdf/1606.08740.pdf</ref> и <ref>https://smartech.gatech.edu/bitstream/handle/1853/55635/PRATAPA-DISSERTATION-2016.pdf</ref>.
+Для исследования масштабируемости была использована параллельная реализация с использованием объекта KSP(решатель СЛАУ) бибилиотеки PETSc(содержащей структуры, типы данных и подпрограммы для параллельной реализации методов на подпространствах Крылова), найденная на GitHub <ref>https://github.com/erdc-cm/petsc-dev/tree/master/src/ksp/ksp/impls/gmres</ref>.
-=== Локальность данных и вычислений ===
+Вычисления проводились на вычислительном кластере, состоящем из 16 узлов, со следующей конфигурацией:
+Altus 1804i Server - 4P Interlagos Node, Quad AMD Opteron 6276, 16C,
+.3 GHz, 128GB, DDR3-1333 ECC, 80GB SSD, MLC, 2.5" HCA, Mellanox ConnectX 2, 1-port QSFP, QDR, memfree, CentOS, Version 5.
-=== Возможные способы и особенности параллельной реализации алгоритма ===
+На рис. 2 представлен график сильной масштабируемости.
-=== Масштабируемость алгоритма и его реализации ===
+[[File:scaling.JPG|thumb|center|500px|'''Рис. 2''' Зависимость производительности от числа процессоров (запуск производился на 1, 8, 27, 64, 216, 512, 1000 процессорах)]]
-=== Динамические характеристики и эффективность реализации алгоритма ===
+== Динамические характеристики и эффективность реализации алгоритма ==
-=== Выводы для классов архитектур ===
+== Выводы для классов архитектур ==
-=== Существующие реализации алгоритма ===
+== Существующие реализации алгоритма ==
 Реализации алгоритма GMRES присутствуют во многих программных пакетах:
@@ Строка 217: / Строка 242: @@
 * [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]
-== Литература ==
+= Литература =
 <references \>

Участник:Alexander34396/Обобщенный метод минимальных невязок: различия между версиями