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

Участник:Анюшева Ирина/Итерационный метод решения системы линейных алгебраических уравнений GMRES (обобщенный метод минимальных невязок): различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показаны 83 промежуточные версии 3 участников)
Строка 1: Строка 1:
−
Основные авторы описания: Анюшева Ирина (614 группа), Исхаков Эльдар (614 группа)
+
{{algorithm
 +
| name              = GMRES (обобщенный метод минимальных невязок)
 +
| serial_complexity = <math>O(m^2n^2)</math>
 +
| pf_height        = <math>O(mn)</math>
 +
| pf_width          = <math>O(mn)</math>
 +
| input_data        = <math>n^2+2n</math>
 +
| output_data      = <math>n</math>
 +
}}
  
−
== Свойства и структура алгоритма ==
+
Основные авторы описания: [[U:Анюшева Ирина|Анюшева Ирина]] (614 группа, разделы: 1.1-1.2, 1.7-1.8, 2.7, 3), Исхаков Эльдар (614 группа, разделы 1.3-1.6, 1.9-1.10)
−
=== Общее описание алгоритма ===
 
−
Обобщенный метод минимальных невязок (GMRES) для решения системы линейных алгебраических уравнений аппроксимирует решение с помощью вектора в подпространства Крылова с минимальным остатком. Метод GMRES разработан Ю.Саадом и Мартин Х. Шульцом в 1986 году. Метод по праву считается одним из самых эффективных численных методов решения несимметричных систем. Он гарантирует неувеличение нормы невязки в ходе итерационного процесса, который строится основан на построении базиса в соответствующей системе подпространстве Крылова <math>K_j(A, r_0)</math>.
 
−
Затем решение уточняется некоторой добавкой, представленной в виде разложения по этому базису.
 
  
−
Метод обобщенных минимальных невязок популярен из-за наличия ряда преимуществ: он ошибкоустойчив, допускает эффективное распараллеливание, не требует нахождения параметра релаксации, обладает суперлинейной скоростью сходимости. Однако в чистом виде для больших систем метод не используется из-за чрезвычайно больших затрат памяти. Зато широкое распространение получила перезапускаемая версия метода, подоразумевающая перезапуск метода каждые m итераций. Уменьшение размерности пространства может привести к стагнации метода – процессу, при котором с каждой следующей итерацией норма невязки уменьшается незначительно, либо не уменьшается вовсе.
+
= Свойства и структура алгоритма =
 +
== Общее описание алгоритма ==
 +
Обобщенный метод минимальных невязок (GMRES) для решения системы линейных алгебраических уравнений аппроксимирует решение с помощью вектора в подпространстве Крылова с минимальной невязкой. Метод GMRES разработан Ю. Саадом и Мартин Х. Шульцом и впервые предложен в 1986 году.  Метод гарантирует неувеличение нормы невязки в ходе итерационного процесса, который основан на построении базиса в соответствующей системе подпространстве Крылова <math>K_j(A, r_0)</math>.
 +
Затем решение уточняется некоторой добавкой, представленной в виде разложения по этому базису. Метод сходится для любого начального приближения.
  
−
=== Математическое описание алгоритма ===
+
Входные данные алгоритма: система линейных алгебраических уравнений, которую необходимо решить.
 +
 
 +
Параметры алгоритма: <math>m</math> - размерность подпространства Крылова, <math>\varepsilon</math> - параметр останова алгоритма.
 +
 
 +
Выходные данные: решение системы в виде вектора.
 +
 
 +
Краткое описание алгоритма:
 +
* Задать параметры алгоритма и начальное приближение;
 +
* Пока не выполнен критерий останова повторить:
 +
** Построить базис подпространства Крылова;
 +
** Минимизировать невязку.
 +
 
 +
== Математическое описание алгоритма ==
 
Исходные данные: система линейных алгебраических уравнений <math>Ax=b</math>, где:
 
Исходные данные: система линейных алгебраических уравнений <math>Ax=b</math>, где:
  
Строка 17: Строка 36:
 
<math>x^*</math> <math>-</math> точное решение системы.
 
<math>x^*</math> <math>-</math> точное решение системы.
  
−
Вычисляемые данные: <math>x^{(s)} - </math>приближенное решение системы.
+
Вычисляемые данные: <math>x^{m} - </math>приближенное решение системы.
  
−
====Ортогонализация Арнольди====
+
В качестве нормы рассматривается евклидова норма.
 +
 
 +
===Ортогонализация Арнольди===
 
Для построения базиса <math>V</math> в подпространстве Крылова в методе GMRES применяется ортогонализация Арнольди.
 
Для построения базиса <math>V</math> в подпространстве Крылова в методе GMRES применяется ортогонализация Арнольди.
  
Строка 32: Строка 53:
 
:<math>
 
:<math>
 
\begin{align}
 
\begin{align}
−
AV_i &= V_i H_i + h_{i+1,i}v_{i+1} e_i^T &=V_{i+1}\overline{H_i},
+
AV_i &= V_i H_i + h_{i+1,i}v_{i+1} e_i^T &=V_{i+1}\overline{H}_i,
 
\end{align}
 
\end{align}
 
</math>
 
</math>
 
:<math>
 
:<math>
 
\begin{align}
 
\begin{align}
−
V_i &= [v_1|v_2|...|v_i], \overline{H_i} - \end{align}
+
V_i &= [v_1|v_2|...|v_i], \overline{H}_i - \end{align}
 
</math> верхняя матрица Хессенберга размерности <math>(i+1)\times i</math>
 
</math> верхняя матрица Хессенберга размерности <math>(i+1)\times i</math>
  
−
====Минимизация невязки====
+
===Минимизация невязки===
  
 
Невязку уравнения можно определить как
 
Невязку уравнения можно определить как
−
<math> J(y) = \parallel b- A x\parallel = \parallel b- (Ax_0+V_m y)\parallel =\parallel \beta e_1 - \overline{H_m} y \parallel</math>,  
+
<math> J(y) = \parallel b- A x\parallel = \parallel b- (Ax_0+V_m y)\parallel =\parallel \beta e_1 - \overline{H}_m y \parallel</math>,  
 
где <math>y -</math> вектор размерности <math>m</math>.
 
где <math>y -</math> вектор размерности <math>m</math>.
  
Строка 49: Строка 70:
 
<math>x_m = x_0 + V_m y_m</math>, вектор <math>y_m</math> может найден как решение линейной задачи наименьших квадратов размера <math>(m+1)\times m</math>, где <math>m<<n</math>.
 
<math>x_m = x_0 + V_m y_m</math>, вектор <math>y_m</math> может найден как решение линейной задачи наименьших квадратов размера <math>(m+1)\times m</math>, где <math>m<<n</math>.
  
−
===Схема реализации последовательного алгоритма===
 
  
−
====Шаги алгоритма====
+
Для решения задачи минимизации приведем матрицу <math>\overline{H}_m</math> к верхнему треугольному виду с помощью вращений Гивенса. На шаге
 +
<math>i</math> матрица вращения <math>\Omega_i</math> имеет размеры <math>(m+1)\times (m+1)</math> и имеет следующий вид:
 +
 
 +
<math>
 +
\Omega_i= \begin{pmatrix}
 +
1 & &&&&&& \\
 +
&\ddots &&&&&&\\
 +
&&1&&&&& \\       
 +
&&&c_i&s_i &&& \\
 +
&&&-s_i& c_i &&&\\
 +
&&&&&1&&& \\
 +
&&&&&&\ddots&\\
 +
&&&&&&&1 \\ 
 +
\end{pmatrix}</math>
 +
 
 +
==Вычислительное ядро алгоритма==
 +
 
 +
В алгоритме GMRES присутствуют два ресурсоемких с точки зрения вычислительных мощностей этапа:
 +
 
 +
* умножение исходной матрицы на вектор при вычислении невязок;
 +
 
 +
* вычисление скалярных произведений векторов при построении ортонормированного базиса подпространства Крылова.
 +
 
 +
==Макроструктура алгоритма==
 +
 
 +
Нижеперечисленные шаги алгоритма и входящие в них операции являются наиболее ресурсоемкими. Они выполняются в приведенном порядке. Каждый из них допускает распараллеливание:
 +
 
 +
* Вычисление нормы невязки для задачи "начальных" условий;
 +
 
 +
* Вычисление элементов Хессенберговой матрицы - матрично-векторные умножения и скалярные произведения векторов;
 +
 
 +
* Вычисление произведений разреженных матриц при минимизации невязки.
 +
 
 +
==Схема реализации последовательного алгоритма==
 +
 
 +
===Шаги алгоритма===
  
 
1. Выбрать произвольное <math>x_0</math>; вычислить <math>r_0=b-Ax_0</math> и <math>v_1=\frac{r_0}{\parallel r_0 \parallel}</math>
 
1. Выбрать произвольное <math>x_0</math>; вычислить <math>r_0=b-Ax_0</math> и <math>v_1=\frac{r_0}{\parallel r_0 \parallel}</math>
Строка 59: Строка 114:
 
3. Для <math>j=1,\dots,m</math> выполнять:
 
3. Для <math>j=1,\dots,m</math> выполнять:
  
−
<math>\bullet \ </math> для <math>i=1,\dots,j</math> выполнять: <math>h_{i,j}=(Av_j,v_i)</math>
+
* для <math>i=1,\dots,j</math> выполнять: <math>h_{i,j}=(Av_j,v_i)</math>
−
<math>\bullet \ \omega_j=Av_j-\sum_{i=1}^j h_{i,j}v_i</math>
+
* <math>\omega_j=Av_j-\sum_{i=1}^j h_{i,j}v_i</math>
−
<math>\bullet \ h_{j+1,j}={\parallel \omega_j \parallel}</math>
+
* <math>h_{j+1,j}={\parallel \omega_j \parallel}</math>
−
<math>\bullet \ v_{j+1}=\frac{\omega_j}{h_{j+1,j}}</math>
+
* <math>v_{j+1}=\frac{\omega_j}{h_{j+1,j}}</math>
  
 
4. Вычислить <math>y_m=\arg\min_{y} \parallel \beta e_1 - \overline{H}_my \parallel,\ y \in \mathbb{R}^m</math> и <math>x_m=x_0+V_my_m</math>
 
4. Вычислить <math>y_m=\arg\min_{y} \parallel \beta e_1 - \overline{H}_my \parallel,\ y \in \mathbb{R}^m</math> и <math>x_m=x_0+V_my_m</math>
  
−
5. Вычислить <math>r_m=f-Ax_m</math>. Если достигнутая точность удовлетворительна, остановиться.
+
5. Вычислить <math>r_m=b-Ax_m</math>. Если достигнутая точность удовлетворительна, остановиться.
  
 
Иначе принять <math>x_0=x_m</math> и <math>v_1=\frac{r_m}{\parallel r_m \parallel}</math> и вернуться к пункту 2.
 
Иначе принять <math>x_0=x_m</math> и <math>v_1=\frac{r_m}{\parallel r_m \parallel}</math> и вернуться к пункту 2.
  
−
====Минимизация невязки====
+
===Решение задачи минимизации невязки===
 +
 
 +
Для нахождения <math>y_m=\arg\min_{y} \parallel \beta e_1 - \overline{H}_my \parallel</math> необходимо выполнить следующие шаги:
 +
 
 +
1. Найти <math>m</math> матриц вращений Гивенса <math>\Omega_i(c_i,s_i),\ i=1,\dots,m</math> размера <math>(m+1)\times(m+1)</math>, где <math>s_i=\frac{h_{i+1,i}}{\sqrt{h_{i,i}^2+h_{i+1,i}^2}}</math>, <math>c_i=\frac{h_{i,i}}{\sqrt{h_{i,i}^2+h_{i+1,i}^2}}</math>
 +
 
 +
2. Вычислить <math>\overline{H}_m^{(m)}=\Omega_1\dots\Omega_m\overline{H}_m</math> и <math>\overline{g}_m=\Omega_1\dots\Omega_m\beta e_1</math>, удалить последнюю строку и последний элемент соответственно и решить СЛАУ с полученной верхнетреугольной матрицей и правой частью. Полученное решение <math>y_m</math> минимизирует исходный функционал.
 +
 
 +
==Последовательная сложность алгоритма==
 +
 
 +
Для выполнения одной итерации алгоритма необходимо произвести:
 +
 
 +
* <math>3m+1</math> делений;
 +
 
 +
* <math>2m+1</math> вычислений квадратного корня;
 +
 
 +
* <math>n^2+\frac{3m+5}{2}n</math> умножений;
 +
 
 +
* <math>O(\frac{m^2n^2}{2})</math> сложений (вычитаний).
 +
 
 +
==Информационный граф==
 +
Метод состоит из двух блоков: ортогонализация и минимизация невязки. При ортогонализации основными операциями являются умножение матрицы на вектор и скалярное умножение векторов, для которых возможен параллелизм, что показано в блоке на темно-синем фоне. В блоке <math>min</math> решается задача минимизации невязки, которая, строго говоря, является отдельной самостоятельной задачей.
 +
[[file:gmres.jpg|thumb|center|700px|Рис.1. Информационная структура алгоритма GMRES]]
 +
 
 +
==Ресурс параллелизма алгоритма==
 +
Для нахождения матрично-векторных произведений при вычислении невязок матрица СЛАУ должна быть подвергнута декомпозиции на строчные блоки, размеры которых определяются требованиями равномерности загрузки. Соответствующим образом подвергается декомпозиции и вектор правой части. Для наиболее эффективного вычисления скалярных произведений и линейных комбинаций векторов базиса каждый из векторов <math>v_j</math> должен быть разбит на блоки равного размера по числу процессов; каждый из процессов при этом хранит соответствующие блоки всех векторов. Обработка вектора решения дублируется на всех процессах.
 +
 
 +
==Входные и выходные данные алгоритма==
 +
 
 +
'''Входные данные:'''
 +
* матрица <math>A</math> размера <math>n\times n</math>;
 +
* вектор-столбец <math>b</math> размера <math>n</math>;
 +
* вектор-столбец <math>x_0</math> размера <math>n</math> - начальное приближение; 
 +
* натуральное число <math>m<<n</math>, определяющее размерность подпространства Крылова.
 +
 
 +
'''Выходные данные:'''
 +
* вектор-столбец <math>x</math> размера <math>n</math> - решение системы.
 +
 
 +
==Свойства алгоритма==
 +
 
 +
* Вычислительная мощность: весьма маленькая и равна <math>m^2</math>, так как <math>m<<n</math>.
 +
 
 +
* Устойчивость: GMRES в этом плане является очень хорошим методом. Теоретически можно получить точное решение не более чем за <math>n</math> итераций. Однако, очень приемлемые результаты получаются и за <math>m<<n</math> итераций. Однако, есть примеры, в которых видно, что первые <math>m-1</math> итераций дают невязку, равную некой константе, а на итерации под номером <math>m</math> получаем точное решение.
 +
 
 +
=Программная реализация алгоритма=
 +
 
 +
==Особенности реализации последовательного алгоритма==
 +
 
 +
==Локальность данных и вычислений==
 +
 
 +
==Возможные способы и особенности параллельной реализации алгоритма==
 +
 
 +
==Масштабируемость алгоритма и его реализации==
 +
Исследование масштабируемости параллельной реализации итерационного метода решения системы линейных алгебраических уравнений GMRES проводилось в тестовой среде суперкомпьютера "Ломоносов". Алгоритм реализован на языке С с использованием библиотеки IntelMPI версии 4.1.0. Сборка производилась при помощи компилятора, разработанным компанией Intel версии 13.1.0.
 +
Набор и границы значений изменяемых параметров запуска:
 +
*число процессов [1 : 128] со значениями, равными степеням двойки;
 +
*размер матрицы [128 : 2048] с шагом 128.
 +
Использованные команды для компиляции и запуска:
 +
 
 +
    mpicc main3.c -o prog
 +
 
 +
    #!/bin/bash
 +
    for a in `seq 0 7`;
 +
    do
 +
    let b=2**$a
 +
    sbatch -p test -n $b impi ./prog
 +
    sleep 600s
 +
    done
 +
 
 +
 
 +
В результате проведенных экспериментов был получен следующий диапазон эффективности реализации алгоритма:
 +
*минимальная эффективность реализации 0.002%;
 +
*максимальная эффективность реализации 4.92%.
 +
На следующих рисунках приведены графики производительности и эффективности выбранной реализации в зависимости от изменяемых параметров запуска. Размерность подпространства Крылова во всех случаях равна 10.
 +
[[File:Anyusheva iskhakovPerformance.jpg|thumb|center|1500px|Рис.2 Изменение производительности выполнения в зависимости от числа процессов и размера входной матрицы]]
 +
[[File:Anyusheva iskhakovEfficiency.jpg|thumb|center|1500px|Рис.3 Изменение эффективности выполнения в зависимости от числа процессов и размера входной матрицы]]
 +
Из рисунков видно, что при увеличении количества процессов, а также при увеличении размерности задачи, производительность параллельной реализации алгоритма уменьшается, что говорит о плохой масштабируемости.
 +
Исследуемая реализация алгоритма доступна по [http://pastebin.com/vEyb3kuS ссылке].
 +
 
 +
==Динамические характеристики и эффективность реализации алгоритма==
 +
 
 +
==Выводы для классов архитектур==
 +
 
 +
==Существующие реализации алгоритма==
 +
* Реализация алгоритма в библиотеке SciPy [https://docs.scipy.org/doc/scipy-0.17.0/reference/generated/scipy.sparse.linalg.gmres.html]
 +
* Реализация алгоритма в библиотеке IML++ [http://math.nist.gov/iml++/gmres.h.txt]
 +
* Библиотека для MATLAB [http://www.mathworks.com/help/matlab/ref/gmres.html]
 +
 
 +
= Литература=
 +
1. Y. Saad, M.H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856–869
 +
 
 +
2. Е.Е. Тыртышников. Методы численного анализа. – Москва: Академия, 2007
  
−
===Задача минимизации невязки===
+
3. М.Ю. Баландин, Э.П. Шурина. Некоторые оценки эффективности параллельных алгоритмов решения СЛАУ на подпространствах Крылова. - Вычислительные технологии, том 3, №1, 1998

Текущая версия на 01:13, 12 декабря 2016


GMRES (обобщенный метод минимальных невязок)
Последовательный алгоритм
Последовательная сложность [math]O(m^2n^2)[/math]
Объём входных данных [math]n^2+2n[/math]
Объём выходных данных [math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы [math]O(mn)[/math]
Ширина ярусно-параллельной формы [math]O(mn)[/math]


Основные авторы описания: Анюшева Ирина (614 группа, разделы: 1.1-1.2, 1.7-1.8, 2.7, 3), Исхаков Эльдар (614 группа, разделы 1.3-1.6, 1.9-1.10)

Содержание

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

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

Обобщенный метод минимальных невязок (GMRES) для решения системы линейных алгебраических уравнений аппроксимирует решение с помощью вектора в подпространстве Крылова с минимальной невязкой. Метод GMRES разработан Ю. Саадом и Мартин Х. Шульцом и впервые предложен в 1986 году. Метод гарантирует неувеличение нормы невязки в ходе итерационного процесса, который основан на построении базиса в соответствующей системе подпространстве Крылова [math]K_j(A, r_0)[/math]. Затем решение уточняется некоторой добавкой, представленной в виде разложения по этому базису. Метод сходится для любого начального приближения.

Входные данные алгоритма: система линейных алгебраических уравнений, которую необходимо решить.

Параметры алгоритма: [math]m[/math] - размерность подпространства Крылова, [math]\varepsilon[/math] - параметр останова алгоритма.

Выходные данные: решение системы в виде вектора.

Краткое описание алгоритма:

  • Задать параметры алгоритма и начальное приближение;
  • Пока не выполнен критерий останова повторить:
    • Построить базис подпространства Крылова;
    • Минимизировать невязку.

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

Исходные данные: система линейных алгебраических уравнений [math]Ax=b[/math], где:

[math]A=(a_{ij})[/math] [math]-[/math] вещественная матрица размера [math]n \times n[/math];

[math]b[/math] и [math]x[/math] вектора из [math]n[/math] элементов;

[math]x^*[/math] [math]-[/math] точное решение системы.

Вычисляемые данные: [math]x^{m} - [/math]приближенное решение системы.

В качестве нормы рассматривается евклидова норма.

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

Для построения базиса [math]V[/math] в подпространстве Крылова в методе GMRES применяется ортогонализация Арнольди.

Формулы метода:

[math] \begin{align} (h_{i+1,i}v_{i+1}) & = Av_i - \sum_{j = 1}^{i} h_{j,i}v_j, \\ v_1 &= \frac{r_0}{\parallel r_0 \parallel}, \end{align} [/math]

В матрично-векторных обозначениях соотношение может быть записано как

[math] \begin{align} AV_i &= V_i H_i + h_{i+1,i}v_{i+1} e_i^T &=V_{i+1}\overline{H}_i, \end{align} [/math]
[math] \begin{align} V_i &= [v_1|v_2|...|v_i], \overline{H}_i - \end{align} [/math] верхняя матрица Хессенберга размерности [math](i+1)\times i[/math]

1.2.2 Минимизация невязки

Невязку уравнения можно определить как [math] J(y) = \parallel b- A x\parallel = \parallel b- (Ax_0+V_m y)\parallel =\parallel \beta e_1 - \overline{H}_m y \parallel[/math], где [math]y -[/math] вектор размерности [math]m[/math].

Вектор [math]x_m[/math] может быть получен как [math]x_m = x_0 + V_m y_m[/math], вектор [math]y_m[/math] может найден как решение линейной задачи наименьших квадратов размера [math](m+1)\times m[/math], где [math]m\lt \lt n[/math].


Для решения задачи минимизации приведем матрицу [math]\overline{H}_m[/math] к верхнему треугольному виду с помощью вращений Гивенса. На шаге [math]i[/math] матрица вращения [math]\Omega_i[/math] имеет размеры [math](m+1)\times (m+1)[/math] и имеет следующий вид:

[math] \Omega_i= \begin{pmatrix} 1 & &&&&&& \\ &\ddots &&&&&&\\ &&1&&&&& \\ &&&c_i&s_i &&& \\ &&&-s_i& c_i &&&\\ &&&&&1&&& \\ &&&&&&\ddots&\\ &&&&&&&1 \\ \end{pmatrix}[/math]

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

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

  • умножение исходной матрицы на вектор при вычислении невязок;
  • вычисление скалярных произведений векторов при построении ортонормированного базиса подпространства Крылова.

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

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

  • Вычисление нормы невязки для задачи "начальных" условий;
  • Вычисление элементов Хессенберговой матрицы - матрично-векторные умножения и скалярные произведения векторов;
  • Вычисление произведений разреженных матриц при минимизации невязки.

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

1.5.1 Шаги алгоритма

1. Выбрать произвольное [math]x_0[/math]; вычислить [math]r_0=b-Ax_0[/math] и [math]v_1=\frac{r_0}{\parallel r_0 \parallel}[/math]

2. Задать матрицу [math]\overline{H}_m[/math] размерности [math](m+1)\times m[/math] и всем ее элементам присвоить нулевые значения

3. Для [math]j=1,\dots,m[/math] выполнять:

  • для [math]i=1,\dots,j[/math] выполнять: [math]h_{i,j}=(Av_j,v_i)[/math]
  • [math]\omega_j=Av_j-\sum_{i=1}^j h_{i,j}v_i[/math]
  • [math]h_{j+1,j}={\parallel \omega_j \parallel}[/math]
  • [math]v_{j+1}=\frac{\omega_j}{h_{j+1,j}}[/math]

4. Вычислить [math]y_m=\arg\min_{y} \parallel \beta e_1 - \overline{H}_my \parallel,\ y \in \mathbb{R}^m[/math] и [math]x_m=x_0+V_my_m[/math]

5. Вычислить [math]r_m=b-Ax_m[/math]. Если достигнутая точность удовлетворительна, остановиться.

Иначе принять [math]x_0=x_m[/math] и [math]v_1=\frac{r_m}{\parallel r_m \parallel}[/math] и вернуться к пункту 2.

1.5.2 Решение задачи минимизации невязки

Для нахождения [math]y_m=\arg\min_{y} \parallel \beta e_1 - \overline{H}_my \parallel[/math] необходимо выполнить следующие шаги:

1. Найти [math]m[/math] матриц вращений Гивенса [math]\Omega_i(c_i,s_i),\ i=1,\dots,m[/math] размера [math](m+1)\times(m+1)[/math], где [math]s_i=\frac{h_{i+1,i}}{\sqrt{h_{i,i}^2+h_{i+1,i}^2}}[/math], [math]c_i=\frac{h_{i,i}}{\sqrt{h_{i,i}^2+h_{i+1,i}^2}}[/math]

2. Вычислить [math]\overline{H}_m^{(m)}=\Omega_1\dots\Omega_m\overline{H}_m[/math] и [math]\overline{g}_m=\Omega_1\dots\Omega_m\beta e_1[/math], удалить последнюю строку и последний элемент соответственно и решить СЛАУ с полученной верхнетреугольной матрицей и правой частью. Полученное решение [math]y_m[/math] минимизирует исходный функционал.

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

Для выполнения одной итерации алгоритма необходимо произвести:

  • [math]3m+1[/math] делений;
  • [math]2m+1[/math] вычислений квадратного корня;
  • [math]n^2+\frac{3m+5}{2}n[/math] умножений;
  • [math]O(\frac{m^2n^2}{2})[/math] сложений (вычитаний).

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

Метод состоит из двух блоков: ортогонализация и минимизация невязки. При ортогонализации основными операциями являются умножение матрицы на вектор и скалярное умножение векторов, для которых возможен параллелизм, что показано в блоке на темно-синем фоне. В блоке [math]min[/math] решается задача минимизации невязки, которая, строго говоря, является отдельной самостоятельной задачей.

Рис.1. Информационная структура алгоритма GMRES

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

Для нахождения матрично-векторных произведений при вычислении невязок матрица СЛАУ должна быть подвергнута декомпозиции на строчные блоки, размеры которых определяются требованиями равномерности загрузки. Соответствующим образом подвергается декомпозиции и вектор правой части. Для наиболее эффективного вычисления скалярных произведений и линейных комбинаций векторов базиса каждый из векторов [math]v_j[/math] должен быть разбит на блоки равного размера по числу процессов; каждый из процессов при этом хранит соответствующие блоки всех векторов. Обработка вектора решения дублируется на всех процессах.

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

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

  • матрица [math]A[/math] размера [math]n\times n[/math];
  • вектор-столбец [math]b[/math] размера [math]n[/math];
  • вектор-столбец [math]x_0[/math] размера [math]n[/math] - начальное приближение;
  • натуральное число [math]m\lt \lt n[/math], определяющее размерность подпространства Крылова.

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

  • вектор-столбец [math]x[/math] размера [math]n[/math] - решение системы.

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

  • Вычислительная мощность: весьма маленькая и равна [math]m^2[/math], так как [math]m\lt \lt n[/math].
  • Устойчивость: GMRES в этом плане является очень хорошим методом. Теоретически можно получить точное решение не более чем за [math]n[/math] итераций. Однако, очень приемлемые результаты получаются и за [math]m\lt \lt n[/math] итераций. Однако, есть примеры, в которых видно, что первые [math]m-1[/math] итераций дают невязку, равную некой константе, а на итерации под номером [math]m[/math] получаем точное решение.

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

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

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

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

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

Исследование масштабируемости параллельной реализации итерационного метода решения системы линейных алгебраических уравнений GMRES проводилось в тестовой среде суперкомпьютера "Ломоносов". Алгоритм реализован на языке С с использованием библиотеки IntelMPI версии 4.1.0. Сборка производилась при помощи компилятора, разработанным компанией Intel версии 13.1.0. Набор и границы значений изменяемых параметров запуска:

  • число процессов [1 : 128] со значениями, равными степеням двойки;
  • размер матрицы [128 : 2048] с шагом 128.

Использованные команды для компиляции и запуска: 

   mpicc main3.c -o prog

   #!/bin/bash 
   for a in `seq 0 7`; 
   do 
   let b=2**$a 
   sbatch -p test -n $b impi ./prog
   sleep 600s 
   done


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

  • минимальная эффективность реализации 0.002%;
  • максимальная эффективность реализации 4.92%.

На следующих рисунках приведены графики производительности и эффективности выбранной реализации в зависимости от изменяемых параметров запуска. Размерность подпространства Крылова во всех случаях равна 10.

Рис.2 Изменение производительности выполнения в зависимости от числа процессов и размера входной матрицы
Рис.3 Изменение эффективности выполнения в зависимости от числа процессов и размера входной матрицы

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

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

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

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

  • Реализация алгоритма в библиотеке SciPy [1]
  • Реализация алгоритма в библиотеке IML++ [2]
  • Библиотека для MATLAB [3]

3 Литература

1. Y. Saad, M.H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856–869

2. Е.Е. Тыртышников. Методы численного анализа. – Москва: Академия, 2007

3. М.Ю. Баландин, Э.П. Шурина. Некоторые оценки эффективности параллельных алгоритмов решения СЛАУ на подпространствах Крылова. - Вычислительные технологии, том 3, №1, 1998

Источник — https://algowiki-project.org/w/ru/index.php?title=Участник:Анюшева_Ирина/Итерационный_метод_решения_системы_линейных_алгебраических_уравнений_GMRES_(обобщенный_метод_минимальных_невязок)&oldid=20676