Метод Ньютона для систем нелинейных уравнений

...пока в процессе работы (Игорь Коньшин)... (используются описания студентов и 4 их реализации)

??? существует также чья-то страница (октябрь 2017) без реализаций ::: https://algowiki-project.org/ru/Метод_Ньютона

Основные авторы статьи: К.О.Шохин и А.А.Лебедев (разделы 1, 2.4.1, 2.7)
Авторы отдельных разделов: А.Чернышев и Н.Захаров (раздел 2.4.2), О.Р.Гирняк и Д.А.Васильков (разделы 1.7, 2.4.3), А.В.Арутюнов и А.С.Жилкин (раздел 2.4.4)

Метод Ньютона решения систем нелинейных уравнений
Последовательный алгоритм
Последовательная сложность	[math]O(L \cdot n^3)[/math], если для решения СЛАУ использовать метод Гаусса, алгоритм сойдется за [math]L[/math] итераций, вычисление значения каждой функции системы и их производных в точке имеет оценку сложности [math]O(n)[/math]
Объём входных данных	[math]n[/math] функций от [math]n[/math] переменных, [math]n^2[/math] частных производных этих функций по каждой переменной, [math]n[/math]-мерный вектор(начальное приближение решения), [math]\varepsilon[/math] - требуемая точность, [math]{\rm max\_iter}[/math] - максимальное число итераций
Объём выходных данных	[math]n[/math]-мерный вектор

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

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

Метод Ньютона решения систем нелинейных уравнений является обобщением метода Ньютона решения нелинейных уравнений, который основан на идее линеаризации. Пусть [math] F(x) : {\mathbb R}^1 \to {\mathbb R}^1[/math] - дифференцируемая функция и необходимо решить уравнение [math]F(x) = 0[/math].

Взяв некоторое [math]x_0[/math] в качестве начального приближения решения, мы можем построить линейную аппроксимацию

[math]F(x)[/math] в окрестности [math]x_0 : F(x_0+h) \approx F(x_0)+F^\prime(x_0)h[/math] и решить получающееся линейное уравнение [math]F(x_0 )+F^\prime (x_0 )h =0[/math].

Таким образом получаем итеративный метод : [math] x_{k+1} = x_k - {F^\prime(x_k)}^{-1}F(x_k), \ k = 0,1,\ldots [/math] .

Данный метод был предложен Ньютоном в 1669 году. Более точно, Ньютон оперировал только с полиномами; в выражении для [math]F(x+h)[/math] он отбрасывал члены более высокого порядка по h , чем линейные. Ученик Ньютона Рафсон в 1690 г. предложил общую форму метода (т. е. не предполагалось что [math]F(x)[/math] обязательно полином и использовалось понятие производной), поэтому часто говорят о методе Ньютона—Рафсона. Дальнейшее развитие исследований связано с именами таких известных математиков, как Фурье, Коши и другие. Например, Фурье доказал в 1818 г., что метод сходится квадратично в окрестности корня, а Коши (1829, 1847) предложил многомерное обобщение метода и использовал метод для доказательства существования решения уравнения.

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

Пусть дана система из [math]n[/math] нелинейных уравнений с [math]n[/math] неизвестными.

[math] \left\{\begin{matrix} f_1(x_1, \ldots, x_n) = 0, \\ f_2(x_1, \ldots, x_n) = 0, \\ \vdots \\ f_n(x_1, \ldots, x_n) = 0. \end{matrix}\right. [/math] , где [math]f_i(x_1, \ldots,x_n) : {\mathbb R}^n \to {\mathbb R}, \ i = 1,\ldots,n[/math] - нелинейные функции, определенные и непрерывно дифференцируемые в некоторой области [math]G \subset {\mathbb R}^n[/math].

Запишем ее в векторном виде:

[math]\overline{x} = {(x_1,x_2,\ldots,x_n)}^T, F(x) ={[f_1(x),f_2(x),\ldots,f_n(x)]}^T, F(x)=0[/math]

Требуется найти такой вектор [math]\overline{x^*} = {(x^*_1,x^*_2,\ldots,x^*_n)}^T[/math], который, при подстановке в исходную систему, превращает каждое уравнение в верное числовое равенство.

При таком подходе формула для нахождения решения является естественным обобщением формулы одномерного итеративного метода:

[math]x^{(k+1)} = x^{(k)} - W^{-1}(x^{(k)})\cdot F(x^{(k)}) , k=0,1,2,\ldots[/math], где

[math] W = \begin{pmatrix} \frac{\partial{f_1(x_1)}}{\partial{x_1}} & \cdots & \frac{\partial{f_1(x_n)}}{\partial{x_n}} \\ \vdots & \ddots & \vdots \\ \frac{\partial{f_n(x_1)}}{\partial{x_1}} & \cdots & \frac{\partial{f_n(x_n)}}{\partial{x_n}} \end{pmatrix}[/math] – матрица Якоби.

В рассмотренных предположениях относительно функции [math]F(\cdot)[/math] при выборе начального приближения [math]x^{(0)}[/math] из достаточно малой окрестности решения [math]\overline{x^*}[/math] имеет место сходимость последовательности [math]\{x^{(k)}\}[/math]. При дополнительном предположении [math]F(\cdot) \in C^2[/math] имеет место квадратичная сходимость метода.

В качестве критерия окончания процесса итераций обычно берут условие [math]\left \| x^{(k+1)} - x^{(k)} \right \| \lt \varepsilon[/math], где [math]\varepsilon[/math] - требуемая точность решения.

Основная сложность метода Ньютона заключается в обращении матрицы Якоби. Вводя обозначение [math]\Delta x^{(k)} = x^{(k+1)} - x^{(k)}[/math] получаем СЛАУ для вычисления [math]\Delta x^{(k)}:[/math] [math]\frac{\partial{F(x^{(k)})}}{\partial{x}} = -F(x^{(k)})[/math]

Тогда [math]x^{(k+1)} = x^{(k)}+ \Delta x^{(k)}.[/math]

Часто метод Ньютона модифицируют следующим образом. По ходу вычислений или заранее выбирают возрастающую последовательность чисел [math]n_0=0, n_1,\ldots[/math]

При [math]n_i \le k \lt n_{i+1}[/math] вычисление [math]\Delta x^{(k)}[/math] осуществляют по следующей формуле:

[math]\frac{\partial{F(x^{n_i})}}{\partial{x}} = -F(x^{(k)}).[/math]

Увеличение числа итераций, сопровождающее такую модификацию, компенсируется «дешевизной» одного шага итерации.

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

Основная вычислительная нагрузка алгоритма заключается в решении СЛАУ:

[math]\frac{\partial F(x^{(k)})}{\partial x}\Delta x^{(k)} = -F(x^{(k)}) \qquad (1)[/math]

Для нахождения [math]\Delta x^{(k)}[/math], по которому вычисляется значение вектора [math]\overline{x}[/math] на очередной итерации: [math]x^{(k+1)} = x^{(k)} + \Delta x^{(k)}.[/math]

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

Как уже было сказано, основную часть каждой итерации метода составляет нахождение матрицы Якоби и решение СЛАУ [math](1)[/math] для нахождения [math]\Delta x^{(k)}.[/math]

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

Формулы метода описаны выше. На рис 1. представлена блок-схема алгоритма, в которой:

1. [math]n[/math] - число уравнений в СЛАУ
2. Max - предельное число итераций
3. [math]\varepsilon[/math] - точность вычислений
4. [math]x^{(0)}[/math] - начальное приближение

Рис.1 Блок-схема последовательного алгоритма

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

Пусть вычисление значения первой производной функции [math]f_i(x_1,x_2,\ldots,x_n)[/math] в точке [math]x^{(k)}[/math] имеет оценку сложности [math]D_i[/math], решение СЛАУ [math](1)[/math] имеет оценку сложности [math]S[/math] (в эту оценку также включаются операции по нахождению элементов матрицы Якоби), а оценка сложности сложения векторов – [math]V[/math]. Тогда оценка сложности одной итерации представляется в виде:

[math]\sum^{n}_{i=1} {D_i} + S + V[/math]

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

В качестве примера можно рассмотреть алгоритм Ньютона, использующий для решения СЛАУ метод Гаусса. Рассмотрим решение системы, для которой справедливы следующие предположения:

• Вычисление значения каждой функции в заданной точке требует [math]O(n)[/math] операций.
• Вычисление значения производной каждой функции в заданной точке требует [math]O(n)[/math] операций.

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

[math]O(n^2) + O(n^3) + O(n) = O(n^3)[/math]

Предположим также, что выполнены все условия для квадратичной сходимости метода, и для решения системы потребуется небольшое (<100) количество итераций.

Пусть алгоритм сойдется за [math]L[/math] итераций, тогда итоговая сложность будет равна [math]O(L \cdot n^3)[/math].

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

Рис.2 Информационный граф алгоритма

На рис. 2 изображена структура информационного графа алгоритма. Блок IN представляет собой выбор начального приближения, необходимой точности. Каждый блок IT соответствует одной итерации алгоритма, блок-схема итерации алгоритма представлена на рис. 3. Блок OUT соответствует успешному окончанию алгоритма.

Как видно, алгоритм является строго последовательным, каждая итерация зависит от результата предыдущей.

Рис.3 Блок-схема итерации алгоритма

Информационный граф метода Ньютона может быть также представлен в виде следующего рисунка:

Рис.4 Информационный граф алгоритма.

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

Хотя сам по себе алгоритм и является строго последовательным, в нем все же присутствует возможность ускорить выполнение за счет распараллеливания отдельных этапов. Основной ресурс параллелизма заключен в решении СЛАУ [math](1)[/math]. Применять можно любой из известных алгоритмов, ориентируясь на известные особенности структуры матрицы Якоби исходной системы.

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

Пусть задача решается на [math]p[/math] процессорах, тогда сложность решения СЛАУ имеет оценку [math]O(\frac{n^3}{p})[/math].

Теперь каждому процессору нужно вычислить не [math]n[/math], а [math]\frac{n}{p}[/math] функций, поэтому для вычисления элементов матрицы Якоби и правой части СЛАУ потребуется [math]O(\frac{n^2}{p})[/math] операций.

Таким образом, получаем следующую оценку параллельной сложности одной итерации алгоритма: [math]O(\frac{n^2}{p}) + O(\frac{n^3}{p}) + O(n) = O(\frac{n^3}{p})[/math].

Пусть алгоритм сходится за [math]L[/math] итераций, тогда итоговая сложность будет равна [math]O(L \cdot \frac{n^3}{p})[/math].

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

При определении объема входных и выходных данных алгоритма будем придерживаться следующих предположений:

• Рассматривается 64-битная система, таким образом каждый адрес занимает 64 бита.
• Все числа, с которыми оперирует алгоритм также занимают 64 бита.
• Каждая функция будет представлена своим адресом, таким образом вклад функций в объем входных данных будет равен [math]n[/math] чисел.
• Производная каждой функции также представлена своим адресом, и вклад производных в объем входных данных будет равен [math]n^2[/math] чисел.

Входными данными алгоритма являются:

1. Функции, образующие систему уравнений.
2. Производные функций, образующих систему уравнений, если их можно задать аналитически.
3. Точность, с которой требуется найти решение.
4. Максимальное число итераций, которое может быть проделано.
5. Начальное приближение [math]x^{(0)}[/math].

Объем входных данных алгоритма равен [math]n ^2 + 2\cdot n + 2 [/math] в случае, когда на вход алгоритму передаются частные производные ([math]n^2[/math] адресов частных производных [math]n[/math] функций по [math]n[/math] переменным, [math]n[/math] адресов функций, образующих систему уравнений, [math]n[/math]-мерный вектор начального приближения, требуемая точность и максимальное количество итераций), когда частные производные приближенно вычисляются непосредственно в программе, объем входных данных равен [math]2\cdot n + 2.[/math]

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

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

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

Соотношение последовательной и параллельной сложностей зависит от того, какие алгоритмы применяются для решения СЛАУ [math](1)[/math]. Результат работы алгоритма сильно зависит от выбора начального приближения решения. При удачном выборе алгоритм достаточно быстро сходится к искомому решению. Глобальная сходимость для многих задач отсутствует.

Существует теорема о достаточном условии сходимости метода Ньютона:

Пусть функция [math]F(x)[/math] непрерывно дифференцируема в открытом выпуклом множестве [math]G \subset {\mathbb R}^n[/math]. Предположим, что существуют [math]\overline{x^*} \in {\mathbb R}^n[/math] и [math]r,\beta \gt 0 [/math], такие что [math]N(\overline{x^*},r) \subset G , F(\overline{x^*}) = 0[/math] , и существует [math]W^{-1}(\overline{x^*})[/math], причем [math]\left \| W^{-1}(\overline{x^*})\right \| \le \beta [/math] и [math]W(x) \in Lip_{\gamma} (N(\overline{x^*},r))[/math]. Тогда существует [math]\varepsilon \gt 0[/math] такое, что для всех [math]x^{(0)} \in N(\overline{x^*}, \varepsilon)[/math] последовательность [math]x^{(1)},x^{(2)},\ldots[/math] порождаемая итерационным процессом сходится к [math]\overline{x^*}[/math] и удовлетворяет неравенству [math]\left \|x^{(k+1)} - \overline{x^*}\right \| \le \beta\cdot\gamma\cdot{\left \| x^{(k)}- \overline{x^*}\right \|}^2[/math].

Использовались следующие обозначения:

• [math]N(x,r)[/math] - открытая окрестность радиуса [math]r[/math] с центром в точке [math]x[/math];

• [math]W(x) \in Lip_{\gamma}(N(x,r))[/math] означает, что [math]W(x)[/math] непрерывна по Липшицу с константой [math]\gamma[/math] в области [math]N(x,r)[/math].

Оставаясь в предположениях, описанных при вычислении объема входных данных, имеем объем входных данных равный [math]n^2 + 2\cdot n + 2 = O(n^2)[/math] чисел.

Оценку числа операций возьмем для реализации метода Ньютона с помощью метода Гаусса: [math]O(n^3)[/math]. Тогда вычислительная мощность есть [math]O(n)[/math].

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

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

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

В настоящем разделе проведено исследование масштабируемости различных параллельных реализаций Метод Ньютона для систем нелинейных уравнений согласно методике AlgoWiki. В основном исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета.

2.4.1 Реализация 1

В исследуемой реализации для решения СЛАУ использовался метод Гаусса. Исследование проводилось для систем размером от 1024 до 4096 уравнений, с шагом 512. Исследуемая система имеет следующий вид:

[math] \left\{\begin{matrix} cos(x_1) - 1 = 0, \\ cos(x_2) - 1 = 0, \\ \vdots \\ cos(x_n) - 1 = 0. \end{matrix}\right. \qquad (2) [/math]

Данная система выбрана так как для нее известно решение, а также начальное приближение [math] x_0 = (0.87, 0.87, ..., 0.87) [/math], при котором метод Ньютона успешно сходится.

Количество ядер процессоров, на которых производились вычисления, изменялось от 16 до 128 по степеням 2. Эксперименты проводились на суперкомпьютере Ломоносов в разделе test. Исследование проводилось на узлах, обладающих следующими характеристиками:

• Количество ядер: 8 ядер архитектуры x86
• Количество памяти: 12Гб
• Количество GPU: 0

На рис. 5 показаны результаты измерений времени решения системы [math](2)[/math] в зависимости от размера системы и количества процессоров. Можно утверждать, что увеличение числа процессоров дает выигрыш во времени, однако, из-за особенностей реализации параллельного решения СЛАУ, при большом количестве узлов и большом размере системы будет осуществляться пересылка между узлами больших объемов данных и возникающие при этом задержки нивелируют выигрыш от увеличения числа процессоров. Поэтому масштабируемость исследуемой реализации весьма ограничена.

Рис.5 Время решения системы уравнений в зависимости от числа процессоров и размера задачи

Использовался компилятор g++, входящий в состав gcc version 4.4.7 20120313 (Red Hat 4.4.7-4). При компиляции указывался оптимизационный ключ -O2, использовалась библиотека Intel MPI 5.0.1. Ссылка на программную реализацию метода Ньютона для решения систем нелинейных уравнений: https://github.com/ArtyLebedev/newton

2.4.2 Реализация 2

В качестве модельной задачи рассматривался один из примеров, поставляемых вместе с модулем SNES пакета PETSc.

Тестирование алгоритма проводилось на суперкомпьютере "Ломоносов"Суперкомпьютерного комплекса Московского университета.

• Версия MPI - impi/4.1.0
• Реализация BLAS - mkl/11.2.0
• Запуски проводились в сегменте regulal4
• Модель используемого CPU - Intel Xeon X5570 2.93GHz
• Версия PETSc - 3.7.4
• Флаги компиляции: -fPIC -Wall -Wwrite-strings -Wno-strict-aliasing -Wno-unknown-pragmas -fvisibility=hidden -g3

Кроме того, использовались следующие параметры:

1. Pестарт алгоритма GMRES через 300 итераций
2. Максимальное число итераций для нахождения решения СЛАУ - 1500
3. Максимальное число итераций метода Ньютона на данном этапе решения - 20

Набор начальных параметров:

• Число процессоров [1 2 4 8 16 32 48 64 80 96 112 128];
• Порядок матрицы [82 122 162 202 242].

Рис. 6. Изменение производительности в зависимости от числа процессоров и размера матрицы.

Рис. 7. Изменение ускорения в зависимости от числа процессоров и размера матрицы.

Рис. 8. Изменение эффективности в зависимости от числа процессоров и размера матрицы

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

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

2.4.3 Реализация 3

Для исследования масштабируемости алгоритм был реализован нами самостоятельно на библиотеке CUDA 8.0 для языка C++(компилятор nvcc) для Windows на домашнем компьютере(OS Windows 10 x64, CPU Intel Core i7-2600 3.4 GHz, 8 GB RAM, GPU NVIDIA GTX 550 Ti). GPU поддерживает спецификацию CUDA 2.1 и максимальное число потоков в блоке 1024. Однако физически у нас есть всего 192 CUDA ядра(см спецификации видеокарты). Все потоки должны быть элементами двумерного квадратного блока. Поэтому исследования масштабируемости были проведены на 4, 16, 36, 64, 100, 144, 196, 256, 324, 400 потоках(на большем количестве потоков проводить исследования бессмысленно ввиду ухудшения производительности). Параметры компиляции:

"nvcc.exe" -gencode=arch=compute_20,code=\"sm_20,compute_20\" --use-local-env --cl-version 2015 -ccbin "C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\bin"  -I"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v8.0\include" -I"C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v8.0\include" -G --keep-dir Debug -maxrregcount=0 --machine 32 --compile -cudart static  -g   -DWIN32 -D_DEBUG -D_CONSOLE -D_MBCS -Xcompiler "/EHsc /W3 /nologo /Od /FS /Zi /RTC1 /MDd " -o Debug\kernel.cu.obj "kernel.cu"

Решение СЛАУ было реализовано в 2 этапа:

1) Нахождение обратной матрицы с помощью метода Гаусса-Жордана. Для этого было написано ядро для CUDA:

__global__ void gaussjordan(float *A, float *I, int n, int i)

2) Умножение обратной матрицы на вектор правых частей. Для этого было написано ядро для CUDA:

__global__ void multiply_by_vec(float *mat, float *vec, float *out, const int N)

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

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

Говоря о масштабируемости данного алгоритма, важно отметить следующее:

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

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

Вывод: для оценки масштабируемости необходимо сравнивать одну итерацию алгоритма, и притом, без учета времени вычисления значения функций и их производных. В итоге, приходим к выводу, что целесообразно измерить время выполнения одного решения СЛАУ. Также стоит отметить, что в нашей реализации при хорошем выборе начального приближения алгоритм сходился за 10-20 итераций.

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

Рис. 9. Время решения СЛАУ в зависимости от размерности системы и количества потоков

Как уже было сказано, из-за наличия только 192 физических процессоров, данный алгоритм показывает лучшее время работы при использовании 144 либо 256 потоков. 144 потока (размер блока 12х12) - наиболее близкое количество потоков к 192 физическим процессорам. Хорошее время работы на 256 потоках объясняется тем, что мы используем стандартный размер блока 16х16, а также многие операции деления или умножения на количество потоков или размер блока в этом случае компилятор может заменить на операции побитового сдвига, которые работают быстрее обычных.

2.4.4 Реализация 4

Характеристики реализации алгоритма сильно зависят от выбранного способа нахождения матрицы Якоби и решения СЛАУ. Для примера рассматривается реализация алгоритма с использованием метода Гаусса на функциях вида:

[math]f_i(x) = cos(x_i) - 1[/math].

Для этих функций можно задать точное значение производной в любой точке:

[math]f_i^' (x) = -sin(x_i)[/math].

Для тестирования программы было решено использовать исключительно технологию MPI в реализации Intel (IntelMPI) без дополнительных. Тесты проводились на суперкомпьютере Ломоносов в разделе test. Исследование проводилось на узлах, обладающих следующими характеристиками:

• Количество ядер: 8 ядер архитектуры x86
• Количество памяти: 12Гб

Строка компиляции: mpicxx _scratch/Source.cpp -o _scratch/out_file

Строка запуска: sbatch -nN -p test impi _scratch/out_file

Результаты тестирования представлены на рисунках 10 и 11, где рис. 10 отображает время работы данной реализации, а рис. 11 - ускорение:

Рис.2 Время решения системы уравнений в зависимости от числа процессоров и размера задачи

Рис.3 Ускорение решения системы уравнений в зависимости от числа процессоров и размера задачи

Из рис. 10 и 11 видно, что увеличение числа задействованных в вычислениях процессоров дает выигрыш во времени, однако, из-за необходимости обмениваться данными при решении СЛАУ методом Гаусса, при вовлечении слишком большого числа процессоров, время, требуемое на обмен данными может превысить время непосредственного вычисления. Этот эффект ограничивает масштабируемость программы. Для подтверждения этого тезиса приведён рис. 12, отдельно показывающий время работы задачи с размерностью матрицы 1024:

Рис. 12. Время решения системы из 1024 уравнений в зависимости от количества задействованных процессоров

Как видно из рис. 12 время выполнения программы на 128 процессорах возрастает по сравнению с временем работы на 64 процессорах. Это связано с тем, что на каждый процессор приходится недостаточно индивидуальной загрузки и основное время работы программы тратится на передачу данных между процессорами. Если же увеличение количества задействованных процессоров не приводит к возникновению этого эффекта, то увеличение количества процессоров в 2 раза ведёт к увеличению скорости работы программы также приблизительно в 2 раза, что хорошо видно на рис. 12.

Исходный код программы

#include <iostream> #include <cmath> #include <fstream> #include <vector> #include <cstdlib> #include <limits> #include <mpi.h> #include <stdio.h> #include <assert.h> typedef std::vector<double> Vec; typedef double(*Function)(const size_t&, const Vec&); typedef double(*Deriv)(const size_t&, const size_t&, const Vec&); typedef std::vector<Function> Sys_; typedef std::vector<std::vector<Deriv> > Derivatives; typedef std::vector<std::vector<double> > Matrix; typedef Vec(*LSSolver)(const Matrix&, const Vec&); typedef Vec(*VecSum)(const Vec&, const Vec&); typedef double(*VecDiff)(const Vec&, const Vec&); struct LS { Matrix A; Vec b; }; size_t system_size = 4096; double func(const size_t& i, const Vec& v) { return cos(v[i]) - 1; } double derivative(const size_t& i, const size_t& j, const Vec& v) { if (i == j) { return -sin(v[i]); } return 0.0; } static void printVec(const Vec& v) { std::cout << "size: " << v.size() << std::endl; for (size_t i = 0; i < v.size(); ++i) { std::cout << v[i] << " "; } std::cout << std::endl; } Vec gauss(const Matrix& A, const Vec& b) { MPI_Status status; int nSize, nRowsBloc, nRows, nCols; int Numprocs, MyRank, Root = 0; int irow, jrow, icol, index, ColofPivot, neigh_proc; double *Input_A, *Input_B, *ARecv, *BRecv; double *Output, Pivot; double *X_buffer, *Y_buffer; double *Buffer_Pivot, *Buffer_bksub; double tmp; MPI_Comm_rank(MPI_COMM_WORLD, &MyRank); MPI_Comm_size(MPI_COMM_WORLD, &Numprocs); if (MyRank == 0) { nRows = A.size(); nCols = A[0].size(); nSize = nRows; Input_B = (double *)malloc(nSize * sizeof(double)); for (irow = 0; irow < nSize; irow++) { Input_B[irow] = b[irow]; } Input_A = (double *)malloc(nSize*nSize * sizeof(double)); index = 0; for (irow = 0; irow < nSize; irow++) for (icol = 0; icol < nSize; icol++) Input_A[index++] = A[irow][icol]; } MPI_Bcast(&nRows, 1, MPI_INT, Root, MPI_COMM_WORLD); MPI_Bcast(&nCols, 1, MPI_INT, Root, MPI_COMM_WORLD); /* .... Broad cast the size of the matrix to all ....*/ MPI_Bcast(&nSize, 1, MPI_INT, 0, MPI_COMM_WORLD); nRowsBloc = nSize / Numprocs; /*......Memory of input matrix and vector on each process .....*/ ARecv = (double *)malloc(nRowsBloc * nSize * sizeof(double)); BRecv = (double *)malloc(nRowsBloc * sizeof(double)); /*......Scatter the Input Data to all process ......*/ MPI_Scatter(Input_A, nRowsBloc * nSize, MPI_DOUBLE, ARecv, nRowsBloc * nSize, MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Scatter(Input_B, nRowsBloc, MPI_DOUBLE, BRecv, nRowsBloc, MPI_DOUBLE, 0, MPI_COMM_WORLD); /* ....Memory allocation of ray Buffer_Pivot .....*/ Buffer_Pivot = (double *)malloc((nRowsBloc + 1 + nSize * nRowsBloc) * sizeof(double)); /* Receive data from all processors (i=0 to k-1) above my processor (k).... */ for (neigh_proc = 0; neigh_proc < MyRank; neigh_proc++) { MPI_Recv(Buffer_Pivot, nRowsBloc * nSize + 1 + nRowsBloc, MPI_DOUBLE, neigh_proc, neigh_proc, MPI_COMM_WORLD, &status); for (irow = 0; irow < nRowsBloc; irow++) { /* .... Buffer_Pivot[0] : locate the rank of the processor received */ /* .... Index is used to reduce the matrix to traingular matrix */ /* .... Buffer_Pivot[0] is used to determine the starting value of pivot in each row of the matrix, on each processor */ ColofPivot = ((int)Buffer_Pivot[0]) * nRowsBloc + irow; for (jrow = 0; jrow < nRowsBloc; jrow++) { index = jrow*nSize; tmp = ARecv[index + ColofPivot]; for (icol = ColofPivot; icol < nSize; icol++) { ARecv[index + icol] -= tmp * Buffer_Pivot[irow*nSize + icol + 1 + nRowsBloc]; } BRecv[jrow] -= tmp * Buffer_Pivot[1 + irow]; ARecv[index + ColofPivot] = 0.0; } } } Y_buffer = (double *)malloc(nRowsBloc * sizeof(double)); /* ....Modification of row entries on each processor ...*/ /* ....Division by pivot value and modification ...*/ for (irow = 0; irow < nRowsBloc; irow++) { ColofPivot = MyRank * nRowsBloc + irow; index = irow*nSize; Pivot = ARecv[index + ColofPivot]; assert(Pivot != 0); for (icol = ColofPivot; icol < nSize; icol++) { ARecv[index + icol] = ARecv[index + icol] / Pivot; Buffer_Pivot[index + icol + 1 + nRowsBloc] = ARecv[index + icol]; } Y_buffer[irow] = BRecv[irow] / Pivot; Buffer_Pivot[irow + 1] = Y_buffer[irow]; for (jrow = irow + 1; jrow < nRowsBloc; jrow++) { tmp = ARecv[jrow*nSize + ColofPivot]; for (icol = ColofPivot + 1; icol < nSize; icol++) { ARecv[jrow*nSize + icol] -= tmp * Buffer_Pivot[index + icol + 1 + nRowsBloc]; } BRecv[jrow] -= tmp * Y_buffer[irow]; ARecv[jrow*nSize + irow] = 0; } } /*....Send data to all processors below the current processors */ for (neigh_proc = MyRank + 1; neigh_proc < Numprocs; neigh_proc++) { /* ...... Rank is stored in first location of Buffer_Pivot and this is used in reduction to triangular form ....*/ Buffer_Pivot[0] = (double)MyRank; MPI_Send(Buffer_Pivot, nRowsBloc*nSize + 1 + nRowsBloc, MPI_DOUBLE, neigh_proc, MyRank, MPI_COMM_WORLD); } /*.... Back Substitution starts from here ........*/ /*.... Receive from all higher processors ......*/ Buffer_bksub = (double *)malloc(nRowsBloc * 2 * sizeof(double)); X_buffer = (double *)malloc(nRowsBloc * sizeof(double)); for (neigh_proc = MyRank + 1; neigh_proc < Numprocs; ++neigh_proc) { MPI_Recv(Buffer_bksub, 2 * nRowsBloc, MPI_DOUBLE, neigh_proc, neigh_proc, MPI_COMM_WORLD, &status); for (irow = nRowsBloc - 1; irow >= 0; irow--) { for (icol = nRowsBloc - 1; icol >= 0; icol--) { /* ... Pick up starting Index .....*/ index = (int)Buffer_bksub[icol]; Y_buffer[irow] -= Buffer_bksub[nRowsBloc + icol] * ARecv[irow*nSize + index]; } } } for (irow = nRowsBloc - 1; irow >= 0; irow--) { index = MyRank*nRowsBloc + irow; Buffer_bksub[irow] = (double)index; Buffer_bksub[nRowsBloc + irow] = X_buffer[irow] = Y_buffer[irow]; for (jrow = irow - 1; jrow >= 0; jrow--) Y_buffer[jrow] -= X_buffer[irow] * ARecv[jrow*nSize + index]; } /*.... Send to all lower processes...*/ for (neigh_proc = 0; neigh_proc < MyRank; neigh_proc++) { MPI_Send(Buffer_bksub, 2 * nRowsBloc, MPI_DOUBLE, neigh_proc, MyRank, MPI_COMM_WORLD); } /*.... Gather the result on the processor 0 ....*/ Output = (double *)malloc(nSize * sizeof(double)); MPI_Gather(X_buffer, nRowsBloc, MPI_DOUBLE, Output, nRowsBloc, MPI_DOUBLE, 0, MPI_COMM_WORLD); /* .......Output vector .....*/ Vec res(nSize); if (MyRank == 0) { for (irow = 0; irow < nSize; irow++) res[irow] = Output[irow]; } return res; } Vec sum(const Vec& lhs, const Vec& rhs) { Vec res(lhs.size()); if (lhs.size() != rhs.size()) { return res; } for (size_t i = 0; i < lhs.size(); ++i) { res[i] = lhs[i] + rhs[i]; } return res; } double diff(const Vec& lhs, const Vec& rhs) { double res = 0.; for (size_t i = 0; i < lhs.size(); ++i) { res += lhs[i] - rhs[i]; } return fabs(res); } class Newtone { public: Newtone(int rank) : m_rank(rank) { m_h = std::numeric_limits<double>::epsilon(); } Vec find_solution(const Sys_& sys, const Vec& start, const Derivatives& d, LSSolver solver, VecSum vec_summator, VecDiff vec_differ, const double& eps, const size_t& max_iter) { size_t iter_count = 1; double diff = 0.; Vec sys_val(sys.size(), 0); if (m_rank == 0) { m_jac.reserve(sys.size()); for (size_t i = 0; i < sys.size(); ++i) { Vec v(sys.size()); m_jac.push_back(v); } compute_jacobian(sys, d, start); for (size_t i = 0; i < sys.size(); ++i) { sys_val[i] = -sys[i](i, start); } } MPI_Barrier(MPI_COMM_WORLD); Vec delta = solver(m_jac, sys_val); Vec new_sol(sys.size()), old_sol(sys.size()); if (m_rank == 0) { new_sol = vec_summator(start, delta); old_sol = start; } MPI_Bcast(new_sol.data(), new_sol.size(), MPI_DOUBLE, 0, MPI_COMM_WORLD); MPI_Bcast(old_sol.data(), old_sol.size(), MPI_DOUBLE, 0, MPI_COMM_WORLD); diff = vec_differ(new_sol, old_sol); while (diff > eps && iter_count <= max_iter) { old_sol = new_sol; if (m_rank == 0) { compute_jacobian(sys, d, old_sol); for (size_t i = 0; i < sys.size(); ++i) { sys_val[i] = -sys[i](i, old_sol); } } MPI_Barrier(MPI_COMM_WORLD); delta = solver(m_jac, sys_val); if (m_rank == 0) { new_sol = vec_summator(old_sol, delta); } MPI_Bcast(new_sol.data(), new_sol.size(), MPI_DOUBLE, 0, MPI_COMM_WORLD); iter_count++; diff = vec_differ(new_sol, old_sol); } return new_sol; } double compute_derivative(const size_t& pos, Function func, const size_t& var_num, const Vec& point) { Vec left_point(point), right_point(point); left_point[var_num] -= m_h; right_point[var_num] += m_h; double left = func(pos, left_point), right = func(pos, right_point); return (right - left) / (2 * m_h); } private: void compute_jacobian(const Sys_& sys, const Derivatives& d, const Vec& point) { for (size_t i = 0; i < sys.size(); ++i) { for (size_t j = 0; j < sys.size(); ++j) { double res_val; res_val = d[i][j](i, j, point); m_jac[i][j] = res_val; } } } double m_h; int m_rank; Matrix m_jac; Vec m_right_part; }; int main(int argc, char** argv) { int rank, size; Sys_ s; Derivatives d(system_size); Vec start(system_size, 0.87); for (size_t i = 0; i < system_size; ++i) { s.push_back(&func); } for (size_t i = 0; i < system_size; ++i) { for (size_t j = 0; j < system_size; ++j) d[i].push_back(&derivative); } MPI_Init(&argc, &argv); MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &size); Newtone n(rank); MPI_Barrier(MPI_COMM_WORLD); double time = MPI_Wtime(); Vec sol = n.find_solution(s, start, d, &gauss, &sum, &diff, 0.0001, 100); MPI_Barrier(MPI_COMM_WORLD); time = MPI_Wtime() - time; double max_time; MPI_Reduce(&time, &max_time, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD); if (rank == 0) { std::ofstream myfile; char filename[32]; snprintf(filename, 32, "out_%ld_%d.txt", system_size, size); myfile.open(filename); for (size_t i = 0; i < sol.size(); ++i) { myfile << sol[i] << " "; } myfile << std::endl; myfile << "Time: " << time << " " << max_time << std::endl; myfile.close(); } MPI_Finalize(); return 0; }

myfile << std::endl; myfile << "Time: " << time << " " << max_time << std::endl; myfile.close(); } MPI_Finalize(); return 0; }</source> |}

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

• Последовательные реализации
  • ALIAS C++. Язык реализации - C++. Распространяется бесплатно, исходные коды и примеры использования можно скачать на сайте.
  • Numerical Recipes. Язык реализации - C++. Исходный код можно найти в секции 9.7 книги Numerical recipes(third edition). Бесплатно доступны предыдущие издания для языков C, Fortran77, Fortran90.
  • Sundials. Язык реализации - C, также есть интерфейс для использования в Fortran-программах. Распространяется по лицензии BSD.
  • Numerical Mathematics- NewtonLib. Язык реализации - C, Fortran. Исходные коды доступны в качестве приложения к книге Peter Deuflhards "Newton Methods for Nonlinear Problems -- Affine Invariance and Adaptive Algorithms".
  • PETSc. Язык реализации - C, есть интерфейс для Java и Python. Распространяется по лицензии BSD.

• Параллельные реализации
  • Sundials. Язык реализации - C, также есть интерфейс для использования в Fortran-программах. Распространяется по лицензии BSD.
  • PETSc. Язык реализации - C, есть интерфейс для Java и Python. Распространяется по лицензии BSD.
  • Построение параллельной реализации метода Ньютона для решения задач оптимизации, описываются в работе V. A. Garanzha, A. I. Golikov, Yu. G. Evtushenko, and M. Kh. Nguen "Parallel Implementation of Newton’s Method for Solving Large-Scale Linear Programs" .
  • Построение параллельной реализации метода Ньютона, а также результаты некоторых экспериментов, описываются в работе Renato N. Elias Alvaro L. G. A. Coutinho Marcos A. D. Martins Rubens M. Sydenstricker "PARALLEL INEXACT NEWTON-TYPE METHODS FOR THE SUPG/PSPG SOLUTION OF STEADY INCOMPRESSIBLE 3D NAVIER-STOKES EQUATIONS IN PC CLUSTERS" .
  • Описание использования и результатов экспериментов с параллельной реализацией метода Ньютона в задаче фильтрации вязкой сжимаемой многофазной многокомпонентной смеси в пористой среде описано в работе К.Ю.Богачева "Эффективное решение задачи фильтрации вязкой сжимаемой многофазной многокомпонентной смеси на параллельных ЭВМ".

3 Литература

<references \>