Участник:Sergey.protserov/Метод Якоби решения СЛАУ

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

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

Метод Якоби -- одношаговый стационарный итерационный метод решения СЛАУ вида [math]Ay = f[/math], где [math] A = \left( \begin{array}{ccc} a_{11} & \dots & a_{1m} \\ \vdots & \ddots & \vdots \\ a_{m1} & \dots & a_{mm} \\ \end{array} \right) [/math], [math] f = \left( \begin{array}{c} f_{1} \\ \vdots \\ f_{m} \\ \end{array} \right) [/math], [math] y = \left( \begin{array}{c} y_{1} \\ \vdots \\ y_{m} \\ \end{array} \right) [/math], [math]\det A \ne 0[/math].

Достаточным условием сходимости метода является свойство строгого диагонального преобладания у матрицы [math]A[/math].

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

Каноническая форма одношагового стационарного итерационного метода имеет вид : [math] B\frac{y^{n+1} - y^{n}}{\tau} + Ay^{n} = f, \quad n = 0,\,1,\,\dots\,, [/math]

где [math]B[/math] — невырожденная матрица [math]m \times m[/math], [math]\tau \in \mathbb{R}[/math], [math]y^{0}[/math] — заданное начальное приближение. Решение исходной СЛАУ находится приближённо посредством последовательных итераций. На [math]n[/math]-ом шаге находится [math]y^{n+1}[/math] — очередное приближение для искомого решения [math]y[/math].

В методе Якоби [math]\tau = 1[/math], [math]B = D[/math], где [math]D[/math] — диагональная матрица, элементы которой совпадают с элементами, стоящими на главной диагонали матрицы [math]A[/math].

Выражение для [math]y^{n+1}[/math] через [math]y^{n}[/math]: [math]y^{n+1} = D^{-1}\left(D-A\right)y^{n} + D^{-1}f[/math].

В поэлементной записи:

[math]y^{n+1}_{i} = \frac{1}{a_{ii}}\left(f_{i} - \sum_{j=1,\,j \ne i}^{m}a_{ij}y^{n}_{j}\right),\quad i = 1,\,\dots,\,m[/math].

В качестве условия окончания итерационного процесса можно использовать условие [math]\left\lVert y^{n+1} - y^{n}\right\rVert \le \varepsilon[/math], где [math]\varepsilon[/math] — заданная точность. Кроме того, можно ограничить максимальное число итераций, задав [math]n_{max}[/math]. Для оценки ошибки можно использовать невязку [math]Ay^{n+1} - f[/math].

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

Основное время работы алгоритма приходится на последовательные вычисления векторов [math]y^{n+1}[/math] по формуле, приведённой в предыдущем пункте, при уже вычисленных в начале работы алгоритма матрице [math]D^{-1}A[/math] и векторе [math]D^{-1}f[/math].

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

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

1. вычисление [math]D^{-1}[/math]
2. вычисление [math]D^{-1}A[/math]
3. вычисление [math]D^{-1}f[/math]
4. вычисление [math]y^{n+1} = y^{n} - D^{-1}Ay^{n} + D^{-1}f[/math]

Макрооперации 1-3 выполняются единожды, и в силу того, что матрица [math]D[/math] — диагональная, занимают лишь незначительную часть времени работы алгоритма.

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

Кроме того, если в качестве критерия завершения работы алгоритма используется условие [math]\left\lVert y^{n+1} - y^{n}\right\rVert \le \varepsilon[/math], требуется также вычислять указанную величину. В дальнейшем при описании алгоритма мы будем предполагать, что этот критерий не используется, а используется завершение работы по достижении максимального числа итераций.

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

1. составить диагональную матрицу [math]D[/math]
2. вычислить [math]D^{-1}[/math]
3. вычислить [math]D^{-1}A[/math]
4. вычислить [math]D^{-1}f[/math]
5. выполнять вычисления по формуле [math]y^{n+1} = y^{n} - D^{-1}Ay^{n} + D^{-1}f,\quad n = 0,\,1,\,\dots,\,n_{max}[/math]

При этом на [math]n[/math]-ом шаге итераций необходимо хранить оба вектора [math]y^{n}[/math], [math]y^{n+1}[/math].

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

Шаг 2 предыдущего пункта требует выполнения [math]m[/math] операций деления вещественных чисел, шаг 3 (с учётом того, что матрица [math]D^{-1}[/math] — диагональная) требует выполнения [math]m^{2}[/math] операций умножения вещественных чисел, аналогично шаг 4 требует [math]m[/math] операций умножения, а каждая итерация шага 5 требует [math]m^{2}[/math] умножений и [math]m^{2} + 2m[/math] сложений/вычитаний.

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

Граф пятого шага алгоритма в случае матрицы порядка [math]2[/math]. Здесь [math]B = D^{-1}A[/math], [math]c = D^{-1}f[/math]. Квадратами обозначены данные (элементы матриц и векторов, стоящие в указанных в скобках позициях), кругами обозначены операции, пунктирными линиями обозначены перезаписи данных (переход от [math]y^{n}[/math] к [math]y^{n+1}[/math] в ходе итераций), элементы данных и идущие от них стрелки нарисованы разными цветами для наглядности.

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

Шаг 2 требует один ярус из [math]m[/math] операций деления, шаг 3 требует один ярус из [math]m^{2}[/math] операций умножения, шаг 4 требует один ярус из [math]m[/math] операций умножения, а каждая итерация шага 5 требует по [math]m[/math] ярусов умножений и сложений (в каждом из ярусов — [math]m[/math] операций) для выполнения умножения матрицы на вектор и ещё два яруса по [math]m[/math] сложений/вычитаний, причём итерации выполняются последовательно.

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

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

1. Вещественная [math]m \times m[/math] матрица [math]A[/math], вообще говоря, плотная
2. Вещественный [math]m[/math]-мерный вектор правой части [math]f[/math]
3. Вещественный [math]m[/math]-мерный вектор начального приближения [math]y^{0}[/math]
4. Максимальное число итераций алгоритма [math]n_{max}[/math]

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

1. Вещественный [math]m[/math]-мерный вектор приближённого решения [math]y^{n_{max}}[/math]

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

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

Зависимость времени работы алгоритма (100 итераций) в секундах от порядка входной матрицы и количества MPI-процессов:

	Кол-во процессов
Порядок матрицы	1	2	4	8	16	32	64
2500	4.063 с	2.072 с	1.044 с	0.533 с	0.283 с	0.186 с	0.150 с
5000	16.205 с	8.283 с	4.124 с	2.059 с	1.084 с	0.600 с	0.381 с
10000	64.708 с	33.071 с	16.331 с	8.145 с	4.094 с	2.208 с	1.192 с
20000	295.386 с	132.868 с	66.455 с	33.137 с	16.553 с	8.329 с	4.524 с

Из приведённых данных видна хорошая слабая масштабируемость алгоритма, а так же тот факт, что увеличение числа процессов с 32 до 64 при входной матрице порядка 2500 уже не даёт значительного выигрыша во времени работы, т.е., что сильная масштабируемость с некоторого момента начинает падать. В исследуемой реализации не осуществляется распараллеливание шагов 2 и 4 алгоритма, т.к. ожидалось, что распределение выполнения [math]2m[/math] операций по процессам не принесёт значительного выигрыша из-за издержек на пересылки множества малых порций данных, а эксперимент показал, что во всех рассмотренных случаях временные затраты на выполнение шагов 2 и 4 в совокупности не превышают [math]4 \cdot 10^{-4}[/math] секунд.

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

Автору статьи не известно о существовании хороших реализаций данного алгоритма.