Алговики - Вклад [ru]

Обсуждение участника:Alexander34396

2016-12-19T09:02:03Z

Kpetrov: /* Пункт 2.4 */

--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK) К вашей работе есть несколько замечаний по содержанию и вопросы по пункту с масштабируемостью. Новая серия замечаний стоит с моей подписью.

----

== Пункт 1.2.3 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''вычислить <math> r_m </math> и остановиться,если требуемая точность была достигнута, иначе повторить.'''

Как здесь предлагается оценивать достижение требуемой точности?

== Пункт 1.4 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''Умножение вектора на скаляр;'''

'''Деление вектора на скаляр.'''

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

== Пункт 1.5 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

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

Более распространенное название – версия GMRES с рестартами. Хорошо бы в этом разделе пояснить, почему «на практике используется» именно такая версия.

== Пункт 1.6 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

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

Нужно пояснить почему сложностью расчета <math> y_m </math> можно пренебречь.

== Пункт 1.8 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''...нахождение взвешенной суммы векторов.''' Линейная комбинация векторов?

* Взвешенная сумма векторов - это линейная комбинация векторов

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

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

* этот абзац поправил

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

Почему данный подход является «наиболее рациональным»? Например, n (размер вектора) может быть очень большим и вектор не будет помещаться в память одного процессора.

* этот подход здесь рассматрвиается как один из возможных - и в зависимости от размера входных данных и от вычислительной машины и процессора можно применять разные методы распределения входных данных

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

Требуются ли сборки данных в таком подходе на каждой итерации метода?

* Сбор данных потребуется на каждой итерации метода

== Пункт 1.10 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''Если используется неперезапускаемый метод, GMRES будет сходиться за не более <math> n </math> шагов.'''

Будет просто сходиться с некоторой точностью, или в точной арифметике – это прямой метод за n шагов?

* Хотя GMRES используется как итерационный метод - теоретически он является прямым: реализуемый в точной арифметике, он приводит к точному решению системы порядка <math>n</math> не более, чем в <math>n</math> шагов

'''Существуют примеры, для которых метод застаивается и сходимость есть только на <math>m</math>-м шаге. Для таких систем, любой выбор <math> m </math> меньше, чем <math> n </math> - не сходится.'''

Т.е. за n-1 шагов ошибка может быть сколь угодно большой? Нужно привести пример, ссылку или пояснить – в каком случае метод не сходится менее чем за n шагов?

* ошибка в таком случае будет больше, чем приемлемая и сходимости, соответственно, не будет; примеров на данный случай к сожалению не смог найти

Опечатка - '''включает в себя чрезмерную работы'''

* опечатку поправил

== Пункт 2.4 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

Запускалась ли данная реализация GMRES в PETSc самостоятельно на указанном кластере?

C чем связано то, что авторы работы [3] из списка литературы получили практически точно такой же график масштабируемости, но для предобусловленного метода GMRES, с помощью той же библиотеки PETSc и на том же самом кластере?

Соответствует ли разбиение входных данных по параллельным процессам в реализации PETSc тому, что приведено в пункте 1.8? Описание алгоритма и то, что тестируется должно быть согласовано. Входит ли в замер времени счета описанная в пункте 1.8 сборка вектора решения на управляющем процессоре?

Нужно также указать размерность задачи, описать структуру матрицы в экспериментах, условие на окончание итерационного процесса.

* самостоятельно данная реализация не запускалась. в других статьях по этой теме есть результаты запуска реализации, взятой с GitHub, с использованием библиотеки MPI. так как из всех статей по одной теме далее будет собираться одна большая, содержащая результаты всех статей, то я решил поискать работы по данному методу, в том числе иностранных научных исследователей, где использовались другие библиотеки, в частности PETSc. код, используемый при запуске, так же можно найти на GitHub, как и реализации метода при помощи других библиотек (в статьях других студентов). в указанных в списке используемой литературы статьях (для данной темы нашлись статьи, в которых использовались предобуславливатели, и поэтому здесь указан график для преобусловленного метода) описаны также задачи, на которых можно применять данный метод и структуры матриц, которые соответствуют системам уравнений. точность метода, то есть допустимая ошибка, равна 1E-6. Про время явно в этих статьях не указано, но в других работах (не совсем подходящих для данной статьи), видел указание, что время замерялось полное, то есть с учетом затрат временных ресурсов на сборку векторов, кроме того в библиотеке PETSc можно обращаться к часам явно через команды PetscLogDouble time; и PetscGetTime(&time);CHKERRQ(ierr). насчет размера системы - в статье указано, что решается линейная система, полученная из дискритезации уравнения Гельмгольца (то есть берется трехмерная сетка с определенным шагом и в каждой точке сетки фиксируют значения) - берется суперклетка размером 3х3х3 и периодически повторяется вдоль одной из осей с факторами 1, 2, 3, 8, 16; в суперклетках равномерно от идеальных позиций распределены атомы, для которых как раз решается система уравнений (7a)-(7b), указанная в статье. размер системы здесь - это получается количество атомов в этих суперклетках.
----

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

Здравствуйте. Немного исправил статью, исходя из ваших замечаний.

= 1.2 =
На мой взгляд, нужны определения
span?
l2-норма?

: [[Участник:Alexander34396|Alexander34396]] Дано пояснение обозначения span.Вместо l2-норма теперь используется Евклидова норма.

= 1.2.1 =
Вот я читаю этот раздел как человек со стороны и вижу много формул. Почему они именно такие? В чем суть этого алгоритма?

: [[Участник:Alexander34396|Alexander34396]] Перед алгоритмом добавлено описание идеи алгоритма и пояснения к формулам.

= 1.2.2 =
Новые определения

: [[Участник:Alexander34396|Alexander34396]] Добавлены некоторые пояснения.Добавлена ссылка на статью - вращения Гивенса.

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

: [[Участник:Alexander34396|Alexander34396]] На мой взгляд, формат статей algo-wiki не подразумевает подробных объяснений и выводов формул, поэтому я добавил только некоторые пояснения, которые делают повествование более последовательным и несколько более понятным. При этом для полного понимания описанного метода для человека со стороны всё равно придётся обращаться к сторонним источникам.

Можете или добавить объяснений или возразить мне. В любом случае, пока с прочтением остановился, жду ответа или исправление (а лучше и того, и другого)

= 2.4 =
[[Участник:Dan|Dan]] ([[Обсуждение участника:Dan|обсуждение]]) 19:09, 30 ноября 2016 (MSK) '''принято'''

Обсуждение участника:Alexander34396

2016-12-19T09:01:30Z

Kpetrov: /* Пункт 2.4 */

--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK) К вашей работе есть несколько замечаний по содержанию и вопросы по пункту с масштабируемостью. Новая серия замечаний стоит с моей подписью.

----

== Пункт 1.2.3 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''вычислить <math> r_m </math> и остановиться,если требуемая точность была достигнута, иначе повторить.'''

Как здесь предлагается оценивать достижение требуемой точности?

== Пункт 1.4 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''Умножение вектора на скаляр;'''

'''Деление вектора на скаляр.'''

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

== Пункт 1.5 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

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

Более распространенное название – версия GMRES с рестартами. Хорошо бы в этом разделе пояснить, почему «на практике используется» именно такая версия.

== Пункт 1.6 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

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

Нужно пояснить почему сложностью расчета <math> y_m </math> можно пренебречь.

== Пункт 1.8 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''...нахождение взвешенной суммы векторов.''' Линейная комбинация векторов?

* Взвешенная сумма векторов - это линейная комбинация векторов

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

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

* этот абзац поправил

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

Почему данный подход является «наиболее рациональным»? Например, n (размер вектора) может быть очень большим и вектор не будет помещаться в память одного процессора.

* этот подход здесь рассматрвиается как один из возможных - и в зависимости от размера входных данных и от вычислительной машины и процессора можно применять разные методы распределения входных данных

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

Требуются ли сборки данных в таком подходе на каждой итерации метода?

* Сбор данных потребуется на каждой итерации метода

== Пункт 1.10 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''Если используется неперезапускаемый метод, GMRES будет сходиться за не более <math> n </math> шагов.'''

Будет просто сходиться с некоторой точностью, или в точной арифметике – это прямой метод за n шагов?

* Хотя GMRES используется как итерационный метод - теоретически он является прямым: реализуемый в точной арифметике, он приводит к точному решению системы порядка <math>n</math> не более, чем в <math>n</math> шагов

'''Существуют примеры, для которых метод застаивается и сходимость есть только на <math>m</math>-м шаге. Для таких систем, любой выбор <math> m </math> меньше, чем <math> n </math> - не сходится.'''

Т.е. за n-1 шагов ошибка может быть сколь угодно большой? Нужно привести пример, ссылку или пояснить – в каком случае метод не сходится менее чем за n шагов?

* ошибка в таком случае будет больше, чем приемлемая и сходимости, соответственно, не будет; примеров на данный случай к сожалению не смог найти

Опечатка - '''включает в себя чрезмерную работы'''

* опечатку поправил

== Пункт 2.4 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

Запускалась ли данная реализация GMRES в PETSc самостоятельно на указанном кластере?

C чем связано то, что авторы работы [3] из списка литературы получили практически точно такой же график масштабируемости, но для предобусловленного метода GMRES, с помощью той же библиотеки PETSc и на том же самом кластере?

Соответствует ли разбиение входных данных по параллельным процессам в реализации PETSc тому, что приведено в пункте 1.8? Описание алгоритма и то, что тестируется должно быть согласовано. Входит ли в замер времени счета описанная в пункте 1.8 сборка вектора решения на управляющем процессоре?

Нужно также указать размерность задачи, описать структуру матрицы в экспериментах, условие на окончание итерационного процесса.

* самостоятельно данная реализация не запускалась. в других статьях по этой теме есть результаты запуска реализации, взятой с GitHub, с использованием библиотеки MPI. так как из всех статей по одной теме далее будет собираться одна большая, содержащая результаты всех статей, то я решил поискать работы по данному методу, в том числе иностранных научных исследователей, где использовались другие библиотеки, в частности PETSc. код, используемый при запуске, так же можно найти на GitHub, как и реализации метода при помощи других библиотек (в статьях других студентов). в указанных в списке используемой литературы статьях (для данной темы нашлись статьи, в которых использовались предобуславливатели, и поэтому здесь указан график для преобусловленного метода) описаны также задачи, на которых можно применять данный метод и структуры матриц, которые соответствуют системам уравнений. точность метода, то есть допустимая ошибка, равна 1E-6. Про время явно в этих статьях не указано, но в других работах (не совсем подходящих для данной статьи), видел указание, что время замерялось полное, то есть с учетом затрат временных ресурсов на сборку векторов, кроме того в библиотеке PETSc можно обращаться к часам явно через команды PetscLogDouble time; и PetscGetTime(&time);CHKERRQ(ierr).
насчет размера системы - в статье указано, что решается линейная система, полученная из дискритезации уравнения Гельмгольца (то есть берется трехмерная сетка с определенным шагом и в каждой точке сетки фиксируют значения) - берется суперклетка размером 3х3х3 и периодически повторяется вдоль одной из осей с факторами 1, 2, 3, 8, 16; в суперклетках равномерно от идеальных позиций распределены атомы, для которых как раз решается система уравнений (7a)-(7b), указанная в статье. размер системы здесь - это получается количество атомов в этих суперклетках.
----

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

Здравствуйте. Немного исправил статью, исходя из ваших замечаний.

= 1.2 =
На мой взгляд, нужны определения
span?
l2-норма?

: [[Участник:Alexander34396|Alexander34396]] Дано пояснение обозначения span.Вместо l2-норма теперь используется Евклидова норма.

= 1.2.1 =
Вот я читаю этот раздел как человек со стороны и вижу много формул. Почему они именно такие? В чем суть этого алгоритма?

: [[Участник:Alexander34396|Alexander34396]] Перед алгоритмом добавлено описание идеи алгоритма и пояснения к формулам.

= 1.2.2 =
Новые определения

: [[Участник:Alexander34396|Alexander34396]] Добавлены некоторые пояснения.Добавлена ссылка на статью - вращения Гивенса.

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

: [[Участник:Alexander34396|Alexander34396]] На мой взгляд, формат статей algo-wiki не подразумевает подробных объяснений и выводов формул, поэтому я добавил только некоторые пояснения, которые делают повествование более последовательным и несколько более понятным. При этом для полного понимания описанного метода для человека со стороны всё равно придётся обращаться к сторонним источникам.

Можете или добавить объяснений или возразить мне. В любом случае, пока с прочтением остановился, жду ответа или исправление (а лучше и того, и другого)

= 2.4 =
[[Участник:Dan|Dan]] ([[Обсуждение участника:Dan|обсуждение]]) 19:09, 30 ноября 2016 (MSK) '''принято'''

Обсуждение участника:Alexander34396

2016-12-19T08:09:48Z

Kpetrov: /* Пункт 1.10 */

--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK) К вашей работе есть несколько замечаний по содержанию и вопросы по пункту с масштабируемостью. Новая серия замечаний стоит с моей подписью.

----

== Пункт 1.2.3 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''вычислить <math> r_m </math> и остановиться,если требуемая точность была достигнута, иначе повторить.'''

Как здесь предлагается оценивать достижение требуемой точности?

== Пункт 1.4 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''Умножение вектора на скаляр;'''

'''Деление вектора на скаляр.'''

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

== Пункт 1.5 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

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

Более распространенное название – версия GMRES с рестартами. Хорошо бы в этом разделе пояснить, почему «на практике используется» именно такая версия.

== Пункт 1.6 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

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

Нужно пояснить почему сложностью расчета <math> y_m </math> можно пренебречь.

== Пункт 1.8 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''...нахождение взвешенной суммы векторов.''' Линейная комбинация векторов?

* Взвешенная сумма векторов - это линейная комбинация векторов

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

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

* этот абзац поправил

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

Почему данный подход является «наиболее рациональным»? Например, n (размер вектора) может быть очень большим и вектор не будет помещаться в память одного процессора.

* этот подход здесь рассматрвиается как один из возможных - и в зависимости от размера входных данных и от вычислительной машины и процессора можно применять разные методы распределения входных данных

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

Требуются ли сборки данных в таком подходе на каждой итерации метода?

* Сбор данных потребуется на каждой итерации метода

== Пункт 1.10 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''Если используется неперезапускаемый метод, GMRES будет сходиться за не более <math> n </math> шагов.'''

Будет просто сходиться с некоторой точностью, или в точной арифметике – это прямой метод за n шагов?

* Хотя GMRES используется как итерационный метод - теоретически он является прямым: реализуемый в точной арифметике, он приводит к точному решению системы порядка <math>n</math> не более, чем в <math>n</math> шагов

'''Существуют примеры, для которых метод застаивается и сходимость есть только на <math>m</math>-м шаге. Для таких систем, любой выбор <math> m </math> меньше, чем <math> n </math> - не сходится.'''

Т.е. за n-1 шагов ошибка может быть сколь угодно большой? Нужно привести пример, ссылку или пояснить – в каком случае метод не сходится менее чем за n шагов?

* ошибка в таком случае будет больше, чем приемлемая и сходимости, соответственно, не будет; примеров на данный случай к сожалению не смог найти

Опечатка - '''включает в себя чрезмерную работы'''

* опечатку поправил

== Пункт 2.4 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

Запускалась ли данная реализация GMRES в PETSc самостоятельно на указанном кластере?

C чем связано то, что авторы работы [3] из списка литературы получили практически точно такой же график масштабируемости, но для предобусловленного метода GMRES, с помощью той же библиотеки PETSc и на том же самом кластере?

Соответствует ли разбиение входных данных по параллельным процессам в реализации PETSc тому, что приведено в пункте 1.8? Описание алгоритма и то, что тестируется должно быть согласовано. Входит ли в замер времени счета описанная в пункте 1.8 сборка вектора решения на управляющем процессоре?

Нужно также указать размерность задачи, описать структуру матрицы в экспериментах, условие на окончание итерационного процесса.

----

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

Здравствуйте. Немного исправил статью, исходя из ваших замечаний.

= 1.2 =
На мой взгляд, нужны определения
span?
l2-норма?

: [[Участник:Alexander34396|Alexander34396]] Дано пояснение обозначения span.Вместо l2-норма теперь используется Евклидова норма.

= 1.2.1 =
Вот я читаю этот раздел как человек со стороны и вижу много формул. Почему они именно такие? В чем суть этого алгоритма?

: [[Участник:Alexander34396|Alexander34396]] Перед алгоритмом добавлено описание идеи алгоритма и пояснения к формулам.

= 1.2.2 =
Новые определения

: [[Участник:Alexander34396|Alexander34396]] Добавлены некоторые пояснения.Добавлена ссылка на статью - вращения Гивенса.

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

: [[Участник:Alexander34396|Alexander34396]] На мой взгляд, формат статей algo-wiki не подразумевает подробных объяснений и выводов формул, поэтому я добавил только некоторые пояснения, которые делают повествование более последовательным и несколько более понятным. При этом для полного понимания описанного метода для человека со стороны всё равно придётся обращаться к сторонним источникам.

Можете или добавить объяснений или возразить мне. В любом случае, пока с прочтением остановился, жду ответа или исправление (а лучше и того, и другого)

= 2.4 =
[[Участник:Dan|Dan]] ([[Обсуждение участника:Dan|обсуждение]]) 19:09, 30 ноября 2016 (MSK) '''принято'''

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

2016-12-19T08:09:32Z

Kpetrov: /* Свойства алгоритма */

Основные авторы описания: [[#Участник: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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

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

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> - поправка решения.

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

Для построения ортонормированного базиса <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>.

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

Невязка исходной системы имеет вид:
: <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> последовательных [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|вращений Гивенса]].

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

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

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

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

[[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. Поэтому здесь применяется аппарат параллельных вычислений для более эффективного использования ресурсов памяти вычислительной машины и для повышения эффективности метода в целом.

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

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

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

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

'''Входные данные:'''
* <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>.

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. 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>; выбирая , когда делать перезапуск - это вопрос опыта.

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

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

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

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

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

Масштабируемость рассматривалась в работах <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,
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 представлен график сильной масштабируемости.

[[File:scaling.JPG|thumb|center|500px|'''Рис. 2''' Зависимость производительности от числа процессоров (запуск производился на 1, 8, 27, 64, 216, 512, 1000 процессорах)]]

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

= Литература =
<references \>

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

2016-12-19T07:35:24Z

Kpetrov: /* Ресурс параллелизма алгоритма */

Основные авторы описания: [[#Участник: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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

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

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> - поправка решения.

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

Для построения ортонормированного базиса <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>.

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

Невязка исходной системы имеет вид:
: <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> последовательных [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|вращений Гивенса]].

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

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

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

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

[[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. Поэтому здесь применяется аппарат параллельных вычислений для более эффективного использования ресурсов памяти вычислительной машины и для повышения эффективности метода в целом.

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

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

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

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

'''Входные данные:'''
* <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>.

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. 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>; выбирая , когда делать перезапуск - это вопрос опыта.

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

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

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

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

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

Масштабируемость рассматривалась в работах <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,
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 представлен график сильной масштабируемости.

[[File:scaling.JPG|thumb|center|500px|'''Рис. 2''' Зависимость производительности от числа процессоров (запуск производился на 1, 8, 27, 64, 216, 512, 1000 процессорах)]]

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

= Литература =
<references \>

Обсуждение участника:Alexander34396

2016-12-19T07:35:17Z

Kpetrov: /* Пункт 1.8 */

--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK) К вашей работе есть несколько замечаний по содержанию и вопросы по пункту с масштабируемостью. Новая серия замечаний стоит с моей подписью.

----

== Пункт 1.2.3 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''вычислить <math> r_m </math> и остановиться,если требуемая точность была достигнута, иначе повторить.'''

Как здесь предлагается оценивать достижение требуемой точности?

== Пункт 1.4 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''Умножение вектора на скаляр;'''

'''Деление вектора на скаляр.'''

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

== Пункт 1.5 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

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

Более распространенное название – версия GMRES с рестартами. Хорошо бы в этом разделе пояснить, почему «на практике используется» именно такая версия.

== Пункт 1.6 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

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

Нужно пояснить почему сложностью расчета <math> y_m </math> можно пренебречь.

== Пункт 1.8 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''...нахождение взвешенной суммы векторов.''' Линейная комбинация векторов?

* Взвешенная сумма векторов - это линейная комбинация векторов

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

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

* этот абзац поправил

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

Почему данный подход является «наиболее рациональным»? Например, n (размер вектора) может быть очень большим и вектор не будет помещаться в память одного процессора.

* этот подход здесь рассматрвиается как один из возможных - и в зависимости от размера входных данных и от вычислительной машины и процессора можно применять разные методы распределения входных данных

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

Требуются ли сборки данных в таком подходе на каждой итерации метода?

* Сбор данных потребуется на каждой итерации метода

== Пункт 1.10 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

'''Если используется неперезапускаемый метод, GMRES будет сходиться за не более <math> n </math> шагов.'''

Будет просто сходиться с некоторой точностью, или в точной арифметике – это прямой метод за n шагов?

'''Существуют примеры, для которых метод застаивается и сходимость есть только на <math>m</math>-м шаге. Для таких систем, любой выбор <math> m </math> меньше, чем <math> n </math> - не сходится.'''

Т.е. за n-1 шагов ошибка может быть сколь угодно большой? Нужно привести пример, ссылку или пояснить – в каком случае метод не сходится менее чем за n шагов?

Опечатка - '''включает в себя чрезмерную работы'''

== Пункт 2.4 ==
--[[Участник:Evgeny Mortikov|Evgeny Mortikov]] ([[Обсуждение участника:Evgeny Mortikov|обсуждение]]) 05:42, 2 декабря 2016 (MSK)

Запускалась ли данная реализация GMRES в PETSc самостоятельно на указанном кластере?

C чем связано то, что авторы работы [3] из списка литературы получили практически точно такой же график масштабируемости, но для предобусловленного метода GMRES, с помощью той же библиотеки PETSc и на том же самом кластере?

Соответствует ли разбиение входных данных по параллельным процессам в реализации PETSc тому, что приведено в пункте 1.8? Описание алгоритма и то, что тестируется должно быть согласовано. Входит ли в замер времени счета описанная в пункте 1.8 сборка вектора решения на управляющем процессоре?

Нужно также указать размерность задачи, описать структуру матрицы в экспериментах, условие на окончание итерационного процесса.

----

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

Здравствуйте. Немного исправил статью, исходя из ваших замечаний.

= 1.2 =
На мой взгляд, нужны определения
span?
l2-норма?

: [[Участник:Alexander34396|Alexander34396]] Дано пояснение обозначения span.Вместо l2-норма теперь используется Евклидова норма.

= 1.2.1 =
Вот я читаю этот раздел как человек со стороны и вижу много формул. Почему они именно такие? В чем суть этого алгоритма?

: [[Участник:Alexander34396|Alexander34396]] Перед алгоритмом добавлено описание идеи алгоритма и пояснения к формулам.

= 1.2.2 =
Новые определения

: [[Участник:Alexander34396|Alexander34396]] Добавлены некоторые пояснения.Добавлена ссылка на статью - вращения Гивенса.

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

: [[Участник:Alexander34396|Alexander34396]] На мой взгляд, формат статей algo-wiki не подразумевает подробных объяснений и выводов формул, поэтому я добавил только некоторые пояснения, которые делают повествование более последовательным и несколько более понятным. При этом для полного понимания описанного метода для человека со стороны всё равно придётся обращаться к сторонним источникам.

Можете или добавить объяснений или возразить мне. В любом случае, пока с прочтением остановился, жду ответа или исправление (а лучше и того, и другого)

= 2.4 =
[[Участник:Dan|Dan]] ([[Обсуждение участника:Dan|обсуждение]]) 19:09, 30 ноября 2016 (MSK) '''принято'''

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

2016-11-30T15:34:56Z

Kpetrov: /* Масштабируемость алгоритма и его реализации */

Основные авторы описания: [[#Участник: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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

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

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> - поправка решения.

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

Для построения ортонормированного базиса <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>.

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

Невязка исходной системы имеет вид:
: <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> последовательных [[Метод_Гивенса_(вращений)_QR-разложения_квадратной_матрицы_(вещественный_вариант)|вращений Гивенса]].

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

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

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

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

[[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. Поэтому при увеличении эффективности использования вычислительных ресурсов и ресурсов памяти вычислительной машины при выполнении этих операций повышается эффективность метода в целом. Достигается это применением при выполнении перечисленных операций аппарата параллельных вычислений.

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

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

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

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

'''Входные данные:'''
* <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>.

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. 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>; выбирая , когда делать перезапуск - это вопрос опыта.

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

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

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

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

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

Масштабируемость рассматривалась в работах <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,
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 представлен график сильной масштабируемости.

[[File:scaling.JPG|thumb|center|500px|'''Рис. 2''' Зависимость производительности от числа процессоров (запуск производился на 1, 8, 27, 64, 216, 512, 1000 процессорах)]]

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

= Литература =
<references \>

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

2016-11-08T10:28:41Z

Kpetrov: /* Масштабируемость алгоритма и его реализации */

{{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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

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

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> l_2 </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> l_2 </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> K_m </math> метод использует процесс Арнольди.

'''Исходные данные:'''
* <math> v_1 </math>, такой что <math> \|v_1\|_2 = 1 </math>;
* матрица <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> \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> 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>.

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

Для решения исходной системы 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> y_m </math> может найден как решение линейной задачи наименьших квадратов размера <math> m+1 </math>-на-<math> m </math>:
: <math> y_m = \arg\min_{y}\| \beta e_1 - \overline{H}_my\|_2 </math>.

Для решения задачи матрица <math> \overline{H}_m </math> приводится к верхнему треугольному виду с помощью <math> m </math> последовательных вращений Гивенса.

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

* Умножение матрицы на вектор;
* Вычисление скалярного произведения векторов;
* Вычисление <math> l_2 </math>-нормы вектора;
* Умножение вектора на скаляр;
* Деление вектора на скаляр.

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

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

[[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. Поэтому при увеличении эффективности использования вычислительных ресурсов и ресурсов памяти вычислительной машины при выполнении этих операций повышается эффективность метода в целом. Достигается это применением при выполнении перечисленных операций аппарата параллельных вычислений.

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

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

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

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

'''Входные данные:'''
* <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>.

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. 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>; выбирая , когда делать перезапуск - это вопрос опыта.

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

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

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

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

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

В статье <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>.

[[File:Intel paragon.JPG|thumb|center|800px|'''Рис. 2''' Intel Paragon]]
[[File:Ibm sp2.JPG|thumb|center|800px|'''Рис. 3''' IBM SP2]]
[[File:Cray t3e.JPG|thumb|center|800px|'''Рис. 4''' Cray T3E]]

Теоретические и экспериментальные результаты показывают, что параллельные реализации <math>GMRES(k)</math> являются достаточно масштабируемыми.

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

= Литература =
<references \>

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

2016-10-27T20:09:04Z

Kpetrov: /* Информационный граф */

Основные авторы описания: [[#Участник: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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

= ЧАСТЬ. Свойства и структура алгоритмов =

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> l_2 </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> l_2 </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> K_m </math> метод использует процесс Арнольди.

'''Исходные данные:'''
* <math> v_1 </math>, такой что <math> \|v_1\|_2 = 1 </math>;
* матрица <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> \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> 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>.

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

Для решения исходной системы 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> y_m </math> может найден как решение линейной задачи наименьших квадратов размера <math> m+1 </math>-на-<math> m </math>:
: <math> y_m = \arg\min_{y}\| \beta e_1 - \overline{H}_my\|_2 </math>.

Для решения задачи матрица <math> \overline{H}_m </math> приводится к верхнему треугольному виду с помощью <math> m </math> последовательных вращений Гивенса.

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

* Умножение матрицы на вектор;
* Вычисление скалярного произведения векторов;
* Вычисление <math> l_2 </math>-нормы вектора;
* Умножение вектора на скаляр;
* Деление вектора на скаляр.

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

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

[[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. Поэтому при увеличении эффективности использования вычислительных ресурсов и ресурсов памяти вычислительной машины при выполнении этих операций повышается эффективность метода в целом. Достигается это применением при выполнении перечисленных операций аппарата параллельных вычислений.

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

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

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

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

'''Входные данные:'''
* <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>.

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. 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>; выбирая , когда делать перезапуск - это вопрос опыта.

= ЧАСТЬ. Программная реализация алгоритма =

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

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

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

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

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

= Литература =
<references \>

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

2016-10-27T20:08:00Z

Kpetrov: /* Информационный граф */

Основные авторы описания: [[#Участник: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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

= ЧАСТЬ. Свойства и структура алгоритмов =

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> l_2 </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> l_2 </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> K_m </math> метод использует процесс Арнольди.

'''Исходные данные:'''
* <math> v_1 </math>, такой что <math> \|v_1\|_2 = 1 </math>;
* матрица <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> \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> 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>.

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

Для решения исходной системы 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> y_m </math> может найден как решение линейной задачи наименьших квадратов размера <math> m+1 </math>-на-<math> m </math>:
: <math> y_m = \arg\min_{y}\| \beta e_1 - \overline{H}_my\|_2 </math>.

Для решения задачи матрица <math> \overline{H}_m </math> приводится к верхнему треугольному виду с помощью <math> m </math> последовательных вращений Гивенса.

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

* Умножение матрицы на вектор;
* Вычисление скалярного произведения векторов;
* Вычисление <math> l_2 </math>-нормы вектора;
* Умножение вектора на скаляр;
* Деление вектора на скаляр.

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

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

[[File: Drawing.png|thumb|left|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. Поэтому при увеличении эффективности использования вычислительных ресурсов и ресурсов памяти вычислительной машины при выполнении этих операций повышается эффективность метода в целом. Достигается это применением при выполнении перечисленных операций аппарата параллельных вычислений.

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

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

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

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

'''Входные данные:'''
* <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>.

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. 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>; выбирая , когда делать перезапуск - это вопрос опыта.

= ЧАСТЬ. Программная реализация алгоритма =

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

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

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

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

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

= Литература =
<references \>

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

2016-10-19T06:18:45Z

Kpetrov: /* Свойства алгоритма */

{{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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

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

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> l_2 </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> l_2 </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> K_m </math> метод использует процесс Арнольди.

'''Исходные данные:'''
* <math> v_1 </math>, такой что <math> \|v_1\|_2 = 1 </math>;
* матрица <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> \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> 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>.

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

Для решения исходной системы 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> y_m </math> может найден как решение линейной задачи наименьших квадратов размера <math> m+1 </math>-на-<math> m </math>:
: <math> y_m = \arg\min_{y}\| \beta e_1 - \overline{H}_my\|_2 </math>.

Для решения задачи матрица <math> \overline{H}_m </math> приводится к верхнему треугольному виду с помощью <math> m </math> последовательных вращений Гивенса.

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

* Умножение матрицы на вектор;
* Вычисление скалярного произведения векторов;
* Вычисление <math> l_2 </math>-нормы вектора;
* Умножение вектора на скаляр;
* Деление вектора на скаляр.

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

=== Информационный граф ===
----
[[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>-ой невязки может происходить параллельно
с <math> (i + 1) </math>-ым шагом алгоритма Арнольди.

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

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

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

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

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

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

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

'''Входные данные:'''
* <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>.

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

Обобщенный метод минимальных невязок предназначен для решения несимметричных систем линейных уравнений. 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>; выбирая , когда делать перезапуск - это вопрос опыта.

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

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

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

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

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

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

== Литература ==
<references \>

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

2016-10-18T18:57:16Z

Kpetrov: /* Свойства алгоритма */

{{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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

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

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> l_2 </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> l_2 </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> K_m </math> метод использует процесс Арнольди.

'''Исходные данные:'''
* <math> v_1 </math>, такой что <math> \|v_1\|_2 = 1 </math>;
* матрица <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> \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> 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>.

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

Для решения исходной системы 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> y_m </math> может найден как решение линейной задачи наименьших квадратов размера <math> m+1 </math>-на-<math> m </math>:
: <math> y_m = \arg\min_{y}\| \beta e_1 - \overline{H}_my\|_2 </math>.

Для решения задачи матрица <math> \overline{H}_m </math> приводится к верхнему треугольному виду с помощью <math> m </math> последовательных вращений Гивенса.

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

* Умножение матрицы на вектор;
* Вычисление скалярного произведения векторов;
* Вычисление <math> l_2 </math>-нормы вектора;
* Умножение вектора на скаляр;
* Деление вектора на скаляр.

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

=== Информационный граф ===
----
[[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>-ой невязки может происходить параллельно
с <math> (i + 1) </math>-ым шагом алгоритма Арнольди.

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

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

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

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

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

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

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

'''Входные данные:'''
* <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>.

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

Вычислительная схема требует работы со следующими объектами: разреженной матрицей исходной СЛАУ — <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-схеме присутствует достаточно самостоятельная подзадача — решение переопределенной СЛАУ малой размерности в смысле наименьших квадратов, которая должна рассматриваться отдельно.

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

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

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

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

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

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

== Литература ==
<references \>

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

2016-10-18T18:16:52Z

Kpetrov: /* Ресурс параллелизма алгоритма */

{{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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

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

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <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> l_2 </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> l_2 </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> K_m </math> метод использует процесс Арнольди.

'''Исходные данные:'''
* <math> v_1 </math>, такой что <math> \|v_1\|_2 = 1 </math>;
* матрица <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> \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> 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>.

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

Для решения исходной системы 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> y_m </math> может найден как решение линейной задачи наименьших квадратов размера <math> m+1 </math>-на-<math> m </math>:
: <math> y_m = \arg\min_{y}\| \beta e_1 - \overline{H}_my\|_2 </math>.

Для решения задачи матрица <math> \overline{H}_m </math> приводится к верхнему треугольному виду с помощью <math> m </math> последовательных вращений Гивенса.

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

В общем виде 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> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода 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>.

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

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

* Умножение матрицы на вектор;
* Вычисление скалярного произведения векторов;
* Вычисление <math> l_2 </math>-нормы вектора;
* Умножение вектора на скаляр;
* Деление вектора на скаляр.

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

=== Информационный граф ===
----
[[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>-ой невязки может происходить параллельно
с <math> (i + 1) </math>-ым шагом алгоритма Арнольди.

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

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

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

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

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

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

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

'''Входные данные:'''
* <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>.

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

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

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

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

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

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

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

== Литература ==
<references \>

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

2016-10-18T10:34:34Z

Kpetrov: /* Ресурс параллелизма алгоритма */

{{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]]).

{{algorithm
| name = Обобщенный метод минимальных невязок
| serial_complexity = <math> O(nm^2) </math>
| input_data = <math> (n + 1)^2 </math>
| output_data = <math> n </math>
}}

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

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

'''Обобщённый метод минимальных невязок''' ({{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> на случай систем с несимметричными матрицами.

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

'''Исходные данные:'''
* система линейных алгебраических уравнений вида <math> Ax = b </math>, где <math> A </math> — невырожденная матрица размера <math>n</math>-на-<math> n </math>;
* <math> x_0 </math> - начальное приближение.
'''Вычисляемые данные:'''
* <math> x_k </math> - приближённое решение исходной системы.

Метод GMRES приближает точное решение исходной системы <math> Ax = b </math> вектором <math> x_k </math>, который минимизирует <math> l_2 </math>-норму невязки <math> r_k = Ax_k - b</math> на подпространстве Крылова:
: <math> K_k = K_k(A,b) = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{k-1}b \} </math>.

На каждой итерации метода решение уточняется поправкой, представленной в виде разложения по <math> l_2 </math>-ортонормированному базису пространства <math> K_k </math>:
: <math> x_k = x_0 + z_k </math>, где <math> x_0 </math> - некоторое начальное приближение, <math> z_k \in K_k </math> - поправка решения.

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

Для построения ортонормированного базиса <math> K_k </math> метод использует процесс Арнольди.

'''Исходные данные:'''
* <math> v_1 </math>, такой что <math> \|v_1\|_2 = 1 </math>;
* матрица <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,k </math>.

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

Выполнять для <math> \quad j=1,\ldots,k </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_k = V_ky_k </math>, где <math> y_k \in \mathbb{R}^k </math> - вектор коэффициентов;
# <math> AV_k = V_kH_k + h_{k+1k}v_{k+1} e_k^T = V_{k+1}\overline{H}_k </math>, где <math> V_k = [v_1|v_2|...|v_k] </math>, а <math> \overline{H}_k </math> - верхняя матрица Хессенберга размерности <math> (j+1) </math> на <math> j </math>.

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

Для решения исходной системы GMRES вычисляет приближённое решение <math> x_k </math>, которое минимизизует функционал невязки:
: <math> J(y) = \|b- A x_k\|_2 = \|b - A(x_0 + V_ky)\|_2 = \|\beta e_1 - \overline{H}_ky_k \|_2 </math>, где <math> \beta = \|r_0\|_2 </math>.

Вектор <math> y_k </math> может найден как решение линейной задачи наименьших квадратов размера <math> k+1 </math>-на-<math> k </math>:
: <math> y_k = \arg\min_{y}\| \beta e_1 - \overline{H}_ky\|_2 </math>.

Для решения задачи матрица <math> \overline{H}_k </math> приводится к верхнему треугольному виду с помощью <math> k </math> последовательных вращений Гивенса.

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

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

# найти ортонормированный базис <math> V_k </math> подпространства <math> K_k </math> с помощью ортогонализации Арнольди;
# найти <math> y_k </math>, минимизирующий <math> \|r_k\|_2 </math>;
# вычислить <math> x_k = x_0 + V_ky_k </math>, где <math> x_0 </math> - начальное приближение;
# вычислить <math> r_k </math> и остановиться,если требуемая точность была достигнута, иначе повторить.

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

Вычислительное ядро последовательной версии метода GMRES состоит из двух частей:
* Вычисление ортонормированного базиса <math> K_k </math>;
* Формирование приближенного решения <math> x_k </math>.

Для вычисления ортонормированного базиса <math> K_k </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,k </math>.

Для нахождения приближённого решения <math> x_k </math> метод использует формулу:
: <math> x_k = x_0 + V_ky_k </math>.

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

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

* Умножение матрицы на вектор;
* Вычисление скалярного произведения векторов;
* Вычисление <math> l_2 </math>-нормы вектора;
* Умножение вектора на скаляр;
* Деление вектора на скаляр.

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

На практике используется перезапускаемая версия метода - 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.

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

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

'''1. Сложность вычисления ортонормированного базиса пространства <math> K_m </math>:'''
:a) Для вычисления <math> j </math>-го вектора базиса <math> K_m, j < 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> мультипликативных операций.

=== Информационный граф ===
----
[[File: Drawing.png|300px|caption]]
----
Здесь '''input''' - входные данные алгоритма (квадратная матрица, вектор правой части, начальное приближение, точность решения);
<math> A_i </math> - шаги алгоритма Арнольди, <math> x_i </math> - нахождение невязки на i-той итерации;
'''if''' - вычисление <math> r_k </math> и в случае достижения нужной точности - остановка.

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

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

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

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

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

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

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

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

'''Входные данные:'''
* <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_k </math> - вектор приближённого решения размера <math> n </math>.

'''Объём выходных данных''': <math> n </math>.

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

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

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

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

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

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

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

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

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

Реализации алгоритма GMRES присутствуют во многих программных пакетах:
* [http://www.netlib.org/linpack/ LINPACK]
* [http://math.nist.gov/tnt/download.html/ JAMA/C++]
* [http://math-atlas.sourceforge.net/ ATLAS]
* [http://software.intel.com/mkl‎/ Intel MLK]
* [http://math.nist.gov/iml++/ IML++]
* [http://netlib.org/lapack/ LAPACK]
* [http://osl.iu.edu/research/mtl/ MTL]
* [http://www-users.cs.umn.edu/~saad/software/SPARSKIT/ SPARSKIT]
* [http://researcher.watson.ibm.com/researcher/view_group.php?id=1426 WSMP]

== Литература ==
<references \>

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

2016-10-17T23:21:59Z

Kpetrov: /* Информационный граф */

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

2016-10-17T23:17:11Z

Kpetrov: /* Информационный граф */

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

2016-10-17T22:41:27Z

Kpetrov: /* Информационный граф */

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

2016-10-17T22:31:27Z

Kpetrov: /* Информационный граф */

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

2016-10-17T22:30:58Z

Kpetrov: /* Информационный граф */

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

2016-10-17T22:30:26Z

Kpetrov: /* Последовательная сложность алгоритма */

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

2016-10-17T22:27:59Z

Kpetrov: /* Информационный граф */

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

2016-10-17T20:11:44Z

Kpetrov: /* Общая схема метода */