Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 28 промежуточных версий 5 участников)
Строка 1: Строка 1:
{{Assignment|Dan}}
+
 
 
Основные авторы описания: [[#Участник: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]]).
 
Основные авторы описания: [[#Участник: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]]).
  
Строка 24: Строка 24:
  
 
Метод GMRES приближает точное решение исходной системы <math> Ax = b </math> вектором <math> x_m </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> \operatorname{span} \, \{ a_1, a_2, \ldots, a_k \} </math> - линейная оболочка системы векторов <math> \{ a_1, a_2, \ldots, a_k \} </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> K_m </math>:  
Строка 35: Строка 35:
 
Идея процесса Арнольди состоит в том, что для нахождения вектора <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} </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> 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> 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> 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>.  
  
 
'''Исходные данные:'''  
 
'''Исходные данные:'''  
Строка 56: Строка 58:
 
'''При введении матричных обозначений можно записать:'''  
 
'''При введении матричных обозначений можно записать:'''  
  
# <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>.
  
 
=== Минимизация функционала невязки ===
 
=== Минимизация функционала невязки ===
  
Для решения исходной системы GMRES вычисляет приближённое решение <math> x_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 = 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> 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-разложения_квадратной_матрицы_(вещественный_вариант)|вращений Гивенса]].
  
 
=== Общая схема метода ===
 
=== Общая схема метода ===
Строка 97: Строка 110:
 
* Вычисление скалярного произведения векторов;
 
* Вычисление скалярного произведения векторов;
 
* Вычисление Евклидовой нормы вектора;
 
* Вычисление Евклидовой нормы вектора;
* Умножение вектора на скаляр;
+
* Умножение вектора на скаляр.
* Деление вектора на скаляр.
 
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
Строка 156: Строка 168:
 
В ходе анализа перезапускаемого метода обобщенных минимальных невязок был выявлен ряд подзадач, выполняющихся неоднократно: умножение матрицы на вектор, вычисление скалярного произведения векторов, нахождение нормы вектора, нормирование вектора, нахождение взвешенной суммы векторов. Наиболее трудоемкой среди перечисленных операций является умножение матрицы на вектор.
 
В ходе анализа перезапускаемого метода обобщенных минимальных невязок был выявлен ряд подзадач, выполняющихся неоднократно: умножение матрицы на вектор, вычисление скалярного произведения векторов, нахождение нормы вектора, нормирование вектора, нахождение взвешенной суммы векторов. Наиболее трудоемкой среди перечисленных операций является умножение матрицы на вектор.
  
Выполнение перечисленных операций занимает большую часть времени работы алгоритма GMRES. Поэтому при увеличении эффективности использования вычислительных ресурсов и ресурсов памяти вычислительной машины при выполнении этих операций повышается эффективность метода в целом. Достигается это применением при выполнении перечисленных операций аппарата параллельных вычислений.
+
Выполнение перечисленных операций занимает большую часть времени работы алгоритма GMRES. Поэтому здесь применяется аппарат параллельных вычислений для более эффективного использования ресурсов памяти вычислительной машины и для повышения эффективности метода в целом.
  
 
При выполнении параллельного алгоритма решения систем линейных алгебраических уравнений необходимо распределить входные данные наиболее рациональным способом. Матрица коэффициентов разбивается на непрерывные горизонтальные полосы. Распределение вектора правой части и начального приближения выполнено по принципу дублирования, т.е. все элементы векторов скопированы на все процессоры. Для рассматриваемой задачи решения СЛАУ такой подход распределения допустим, т.к. векторы содержат столько же элементов, сколько и одна строка или один столбец матрицы коэффициентов.
 
При выполнении параллельного алгоритма решения систем линейных алгебраических уравнений необходимо распределить входные данные наиболее рациональным способом. Матрица коэффициентов разбивается на непрерывные горизонтальные полосы. Распределение вектора правой части и начального приближения выполнено по принципу дублирования, т.е. все элементы векторов скопированы на все процессоры. Для рассматриваемой задачи решения СЛАУ такой подход распределения допустим, т.к. векторы содержат столько же элементов, сколько и одна строка или один столбец матрицы коэффициентов.
Строка 162: Строка 174:
 
Подзадачи можно разделить на 2 группы: операции между матрицей и вектором (например, умножение матрицы на вектор) и операции между двумя векторами (например, скалярное умножение). Подзадачи умножения константы на вектор и определение нормы вектора можно относить ко второй группе, так как их трудоемкость одного порядка.
 
Подзадачи можно разделить на 2 группы: операции между матрицей и вектором (например, умножение матрицы на вектор) и операции между двумя векторами (например, скалярное умножение). Подзадачи умножения константы на вектор и определение нормы вектора можно относить ко второй группе, так как их трудоемкость одного порядка.
  
Результатом выполнения рассматриваемых операций является, вектор, части которого расположены на участвовавших в вычислениях процессорах. Для объединения результатов расчета и получения полного вектора на каждом процессоре вычислительной системы необходимо выполнить операцию сбора данных.
+
Результатом выполнения рассматриваемых операций является вектор, части которого расположены на участвовавших в вычислениях процессорах. Для объединения результатов расчета и получения полного вектора на каждом процессоре вычислительной системы необходимо выполнить операцию сбора данных.
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
Строка 189: Строка 201:
 
Если используется неперезапускаемый метод, 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> 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>; выбирая , когда делать перезапуск - это вопрос опыта.
+
Основным недостатком GMRES является то, что объемы работ и ресурсы памяти, требуемые для итерации, линейно возрастают с увеличением числа итераций. Если не удалось получить чрезвычайно быструю сходимость - затраты будут очень высокими. Основной способ преодолеть это ограничение - это перезапуск итерации. После выбранного числа итераций <math> m </math> накопленные данные будут очищены и промежуточные результаты используются в качестве входных данных для последующих итераций. Эта процедура повторяется до тех пор, пока сходимость не будет достигнута. Трудность заключается в выборе подходящего значения <math> m </math>. Если <math> m </math> слишком мало, <math>GMRES(m)</math>  может быть медленно сходиться, или не сойтись совсем. Значение <math> m </math>, большее, чем это необходимо, включает в себя чрезмерную работу и использует больше памяти. К сожалению, нет определенных правил, регулирующих выбор <math> m </math>; выбирая , когда делать перезапуск - это вопрос опыта.
  
 
= Программная реализация алгоритма =
 
= Программная реализация алгоритма =
Строка 201: Строка 213:
 
== Масштабируемость алгоритма и его реализации ==
 
== Масштабируемость алгоритма и его реализации ==
  
В статье <ref> M. Sosonkina, Donald C. S. Allison, Layne T. Watson: Scalability analysis of parallel GMRES implementations </ref> анализируется параллельная реализация <math>GMRES(k)</math> на Intel Paragon, IBM SP2 и Cray T3E.
 
Системы Intel Paragon и IBM SP2 не имеют двухуровневого кэша, что значительно сказывается на времени работы из-за больших потерь скорости и частых обращений к памяти.
 
В Cray T3E каждый процессор соединен с 8-килобайтовым кэшем данных и 96-килобайтовым двухуровневым кэшем данных. Также имеется дополнительный кэш объемом 4 мегабайта. Это достаточная емкость кэш-памяти для хранения всех данных с плавающей точкой с высоким коэффициентом попаданий кэша.
 
 
Ниже представлены графики экспериментальных результатов, где штриховыми линиями указывается теоретическая связь между коэффициентом масштабируемости <math>K</math> и числом процессоров <math>p</math> для фиксированного значения эффективности. Сплошная линия является аппроксимацией экспериментальных данных методом наименьших квадратов. На всех графиках вертикальная ось представляет <math>K</math>, а горизонтальная ось - <math>p</math>.
 
  
 +
Масштабируемость рассматривалась в работах <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>.
  
[[File:Intel paragon.JPG|thumb|center|800px|'''Рис. 2''' Intel Paragon]]
+
Вычисления проводились на вычислительном кластере, состоящем из 16 узлов, со следующей конфигурацией:
[[File:Ibm sp2.JPG|thumb|center|800px|'''Рис. 3''' IBM SP2]]
+
Altus 1804i Server - 4P Interlagos Node, Quad AMD Opteron 6276, 16C,
[[File:Cray t3e.JPG|thumb|center|800px|'''Рис. 4''' Cray T3E]]
+
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 представлен график сильной масштабируемости.
  
Теоретические и экспериментальные результаты показывают, что параллельные реализации <math>GMRES(k)</math> являются достаточно масштабируемыми.
+
[[File:scaling.JPG|thumb|center|500px|'''Рис. 2''' Зависимость производительности от числа процессоров (запуск производился на 1, 8, 27, 64, 216, 512, 1000 процессорах)]]
  
 
== Динамические характеристики и эффективность реализации алгоритма ==
 
== Динамические характеристики и эффективность реализации алгоритма ==

Текущая версия на 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]:

  1. [math] h_{ij} = (Av_j, v_i), \quad i=1,\ldots,j [/math];
  2. [math] \hat{v}_{j+1} = Av_j - \sum_{i=1}^j h_{ij}v_{i} [/math];
  3. [math] h_{j+1j} = \|\hat{v}_{j+1}\|_2 [/math];
  4. [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 может быть записан следующим образом:

  1. найти ортонормированный базис [math] V_m [/math] подпространства [math] K_m [/math] с помощью ортогонализации Арнольди;
  2. найти [math] y_m [/math], минимизирующий [math] \|r_m\|_2 [/math];
  3. вычислить [math] x_m = x_0 + V_my_m [/math], где [math] x_0 [/math] - начальное приближение;
  4. вычислить [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 Литература

<references \>

  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. https://arxiv.org/pdf/1606.08740.pdf
  4. https://smartech.gatech.edu/bitstream/handle/1853/55635/PRATAPA-DISSERTATION-2016.pdf
  5. https://github.com/erdc-cm/petsc-dev/tree/master/src/ksp/ksp/impls/gmres