Участник:Арутюнов А.В.

Эта работа прошла предварительную проверку Дата последней правки страницы: 06.02.2017 Данная работа соответствует формальным критериям. Проверено ASA.

Решение системы нелинейных уравнений методом Ньютона
Последовательный алгоритм
Последовательная сложность	[math]O(n^3) [/math] - одна итерация
Объём входных данных	[math]n^2 + n[/math] функций от [math]n[/math] переменных, [math]n[/math] мерный вектор, [math]x^0[/math] - начальное приближение, ε - точность решения.
Объём выходных данных	[math]n[/math]-мерный вектор вещественных чисел

Автор описания: Арутюнов А.В., Жилкин А.С.

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

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

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

Метод Ньютона для решение систем нелинейных уравнений - обобщение классического метода Ньютона Это итерационный численный метод нахождения корня (нуля) заданной функции. Модификацией метода является метод хорд и касательных.

Классический метод Ньютона или касательных заключается в том, что если [math]x_n[/math] — некоторое приближение к корню [math]x_*[/math] уравнения [math]f(x)=0[/math], [math]f(x) \in C^1[/math], то следующее приближение определяется как корень касательной [math] f(x)[/math] к функции, проведенной в точке [math]x_n[/math].

Метод был предложен Исааком Ньютоном в 1669 году. Поиск решения осуществляется путём построения последовательных приближений и основан на принципах простой итерации.

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

Идея метода заключается в линеаризации уравнений системы

[math] \left\{\begin{matrix}f_1(x_1, ..., x_n) = 0 \\ f_2(x_1,...x_n)=0, \\ ..., \\ f_n(x_1, ..., x_n) = 0, \end{matrix}\right. [/math]

Что позволяет свести исходную задачу СНУ(система нелинейных уравнений) к многократному решению системы линейных уравнений.

Пусть известно приближение [math](x_i)^{(k)}[/math] решения системы нелинейных уравнений [math](x_i)^*[/math].Введем в рассмотрение поправку [math]\Delta x_i[/math] как разницу между решением и его приближением: [math] \Delta x_i = (x_i)^* -(x_i)^{(k)} \Rightarrow (x_i)^*=(x_i)^{(k)} + \Delta x_i, i=\overline{(1,n)} [/math]

Подставим полученное выражение для [math](x_i)^*[/math] в исходную систему.

[math] \left\{\begin{matrix} f_1((x_1)^{(k)} + \Delta x_1, (x_2)^{(k)}+ \Delta x_2, ..., (x_n)^{(k)} + \Delta x_n) = 0, \\ f_2((x_1)^{(k)} + \Delta x_1, (x_2)^{(k)} + \Delta x_2, ..., (x_n)^{(k)} + \Delta x_n) = 0, \\ ... \\ f_n((x_1)^{(k)} + \Delta x_1, (x_2)^{(k)} + \Delta x_2, ..., (x_n)^{(k)} + \Delta x_n) = 0, \end{matrix}\right. [/math]

Неизвестными в этой системе нашейных уравнений являются поправки [math]\Delta x_i [/math]. Для определения [math]\Delta x_i [/math] нужно решить эту систему. Но решить эту задачу так же сложно, как и исходную. Однако эту систему можно линеаризовать и, решив её, получить приближённые значения поправок [math]\Delta x_i[/math] для нашего приближения, т.е. [math]\Delta (x_i)^{(k)}[/math]. Эти поправки не позволяют сразу получить точное решение [math](x_i)^*[/math], но дают возможность приблизиться к решению, - получить новое приближение решения [math]((x_i)^{(k+1)} = ((x_i)^{(k)} + ( \Delta x_i)^{(k)}, i = \overline{(1,n)}[/math]

Для линеаризации системы следует разложить функцию [math]f_i[/math] в ряды Тейлора в окрестности [math](x_i)^k[/math], ограничиваясь первыми дифференциалами. Полученная система имеет вид:

[math]\sum_{i=1}^n\frac{ \partial f_i({x_1}^{(k)}, {x_2}^{(k)}, ..., {x_n}^{(k)})}{ \partial x_i} \Delta {x_i}^{(k)}= -f_j({x_1}^{(k)}, {x_2}^{(k)}, ..., {x_n}^{(k)}), j=\overline{(1,n)}[/math]

Все коэффициенты этого уравнения можно вычислить, используя последнее приближение решения [math](x_i)^{(k)}[/math]. Для решения системы линейных уравнений при [math]n=2,3[/math] можно использовать формулы Крамера, при большей размерности системы [math]n[/math] - метод исключения Гаусса.

Значения поправок используются для оценки достигнутой точности решения. Если максимальная по абсолютной величине поправка меньше заданной точности [math]\varepsilon[/math], расчет завершается. Таким образом, условие окончания расчета:

[math]\delta =\min_{ i=\overline{(1,n)}}|\Delta {x_i}^{(k)}|[/math]

Можно использовать и среднее значение модулей поправок:

[math]\delta = \frac{1}{n}\sum_{i=1}^n |\Delta x_i|\lt \varepsilon [/math]

В матричной форме систему можно записать как:

[math]W(\Delta X^k)*X^{(k)} = -F(X^k)[/math]

где [math]W(x)[/math] - матрица Якоби(производных):

[math]W(x)=(\frac{\partial f_j}{\partial x_i})_{n,n}= \left\{\begin{matrix}(\frac{\partial f_1}{\partial x_1} \frac{\partial f_1}{\partial x_2} ... \frac{\partial f_1}{\partial x_n}) \\ ... ... ... ... \\( \frac{\partial f_n}{\partial x_1} \frac{\partial f_n}{\partial x_2} ... \frac{\partial f_n}{\partial x_n} ) \end{matrix}\right. [/math]

[math]\Delta X^{(k)}= \left\{\begin{matrix} \Delta (x_1)^{(k)} \\ \Delta (x_2)^{(k)} \\ ... \\ \Delta (x_n)^{(k)} \end{matrix}\right. [/math]

[math]F(x)[/math] - вектор-функция

[math]W(X^{(k)})[/math] - матрица Якоби, вычисленная для очередного приближения [math]F(X^{(k)})[/math] - вектор-функция, вычисленная для очередного приближения

Выразим вектор поправок [math]X^{(k)}=-W^{-1}(X^{(k)})*F(X^{(k)})[/math] :

[math]W^{-1}[/math], где [math]W^{-1}[/math]

Окончательная формула последовательных приближений метода Ньютона решения СНУ в матричной форме имеет вид:

[math]X^{(k+1)}=X^{(k)}=W^{-1}(X^{(k)})*F(X^{(k)})[/math]

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

Основная вычислительная нагрузка приходится на

1) Решение СЛАУ: [math]F(X^{(k)})=\frac{\partial F(x^{(k)})}{\partial x} \Delta x^{(k)}[/math]

2)Численное вычисление Якобиана(если производные не даны): [math]\frac{\partial F(x^{(k)})}{\partial x}[/math]

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

Как уже было показано выше, макроструктура алгоритма состоит из двух основных операций в каждой итерации, это нахождением элементов матрицы Якоби и решение СЛАУ.

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

1)Задаётся размерность системы [math]n[/math], требуемая точность, начальное приближённое решение [math]X=(x_i)_n[/math]

2)Вычисляются элементы матрицы Якоби [math]W=\left( {\partial f_i\over \partial x_i} \right)_{n,n}[/math]

3)Вычисляется обратная матрица [math]W^{-1} [/math]

4)Вычисляются вектор функция [math]F=(f_i)_n, f_i=f_i(x_1, x_2,..., x_n), i=1,...,n [/math]

5)Вычисляются вектор поправок [math] \Delta X=W_{-1}*F [/math]

6)Уточняется решение [math]X_{n+1}=X_n+\Delta X[/math]

7)Оценивается достигнутая точность [math]\delta = \max_{i=1,n} \Delta x_i^k[/math]

8)Проверяется условие завершения итерационного процесса [math]\delta\lt =\varepsilon[/math]

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

Сперва стоит отметить, что итоговая сложность алгоритма зависит от того, насколько быстро он будет сходиться. Это в свою очередь зависит от заданного начального приближения [math]x^0[/math] и от условия остановки алгоритма [math]\varepsilon[/math]. Можно однако вычислить сложность одного шага итерации алгоритма.

Предполагаем, что у нас нет значений производных заданных функций

[math]f_1(.), f_2(.), ..., f_n(.)[/math].

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

[math]f_i'(x) = (f(x+h) - f(x-h))/2h [/math],

мы находим приближённое значение производной в интересующей нас точке за 5 операций.

Учитывая, что в Якобиане содержится [math]n^2[/math] элементов - производные каждой функции по каждой переменной, - то для нахождения Якобиана нам суммарно требуется

[math]O(5n^2) = O(n^2)[/math]

операций. Сложность вычислений обратной матрицы

[math]W^{-1} [/math]

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

[math]O(n^3)[/math].

Таким образом сложность вычислительного ядра алгоритма составляет

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

для одного шага алгоритма.

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

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

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

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

Сложность вычисление элементов матрицы Якоби - производных [math]\frac{\partial F(x^{(k)})}{\partial x}[/math] в случае, если они не заданы будет - [math]O(n^2/p)[/math].

Решение СЛАУ [math]F(X^{(k)}) = \frac{\partial F(x^{(k)})}{\partial x} \Delta x^{(k)}[/math] одним из параллельных методов имеет сложность [math]O(n^3/p)[/math].

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

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

Входные данные: [math]n^2 + n[/math] функций от [math]n[/math] переменных, [math]n[/math] мерный вектор, [math]x^0[/math] - начальное приближение, [math]\varepsilon[/math] - точность решения.

Выходные данные: [math]n[/math] вещественных чисел

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

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

Вычислительная мощность алгоритма равна отношению числа операций к суммарному объему входных и выходных данных. Основываясь на описанных выше входных и выходных данных, их общий объём будет равен [math]O(n^2)[/math]. При реализации метода при помощи метода Гаусса число операций будет равно [math]O(n^3)[/math]. Тогда вычислительная мощность алгоритма будет равна [math]O(n)[/math].

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

2.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

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

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

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

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

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

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

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

#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; }

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

Sundials. Язык реализации - C, также есть интерфейс для использования в Fortran-программах. Распространяется по лицензии BSD.

PETSc. Язык реализации - C, есть интерфейс для Java и Python. Распространяется по лицензии BSD.

3 Литература

<references \> https://ru.wikipedia.org/wiki/Метод_Ньютона Тыртышников Е. Е. "Методы численного анализа" — М., Академия, 2007. - 320 c. http://www.intuit.ru/studies/courses/4447/983/lecture/14931

https://software.intel.com/en-us/intel-mpi-library