Метод релаксации матрицы

Метод релаксации матрицы
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(it\_max \cdot N^2) }[/math]
Объём входных данных	[math]\displaystyle{ N^2 }[/math] (матрица значений)
Объём выходных данных	[math]\displaystyle{ N^2 }[/math] (матрица значений)
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ O(it\_max) }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ O(N^2) }[/math]

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

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

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

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

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

Исходная задача: решение уравнения Лапласа [math]\displaystyle{ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0 }[/math] в прямоугольной области [math]\displaystyle{ [0,1]\times[0,1] }[/math] с нулевыми граничными условиями.

Дискретизация на равномерной сетке [math]\displaystyle{ N \times N }[/math] с шагом [math]\displaystyle{ h = 1/(N-1) }[/math] приводит к системе линейных уравнений, которая решается методом релаксации:

Для каждого внутреннего узла [math]\displaystyle{ (i,j) }[/math] на итерации [math]\displaystyle{ k }[/math]: [math]\displaystyle{ u_{i,j}^{(k)} = \frac{u_{i-1,j}^{(k-1)} + u_{i+1,j}^{(k-1)} + u_{i,j-1}^{(k-1)} + u_{i,j+1}^{(k-1)}}{4} }[/math]

Критерий остановки: [math]\displaystyle{ \max_{i,j} |u_{i,j}^{(k)} - u_{i,j}^{(k-1)}| \lt \varepsilon_{\text{max}} }[/math]

Начальное приближение для внутренних точек: [math]\displaystyle{ u_{i,j}^{(0)} = 1 + i + j }[/math]

Контрольная сумма для верификации: [math]\displaystyle{ S = \frac{1}{N^2} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} u_{i,j} \cdot (i+1) \cdot (j+1) }[/math]

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

Вычислительное ядро алгоритма состоит из операций обновления значений во внутренних узлах сетки: [math]\displaystyle{ \text{new_val} = \frac{\text{left} + \text{right} + \text{top} + \text{bottom}}{4} }[/math]

Для каждой внутренней точки на каждой итерации выполняется:

• 4 чтения значений соседей
• 3 операции сложения
• 1 операция деления
• 1 запись нового значения
• 1 операция вычисления погрешности (разности между старым и новым значением)

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

Алгоритм имеет итерационную структуру:

1. Инициализация матрицы (граничные условия + начальное приближение)
2. Для каждой итерации:
   1. Вычисление новых значений для всех внутренних точек (релаксация)
   2. Обновление матрицы
   3. Вычисление максимальной погрешности
   4. Проверка условия сходимости
3. Вычисление контрольной суммы

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

• Инициализация матрицы A граничными условиями и начальными значениями, инициализация матрицы B нулями.
• Цикл по итерациям:

  ```c
  for (it = 1; it <= itmax; it++) {
      eps = 0.0;
      relax();     // Вычисление B[i][j] как среднее соседей из A
      resid();     // Обновление A[i][j] = B[i][j], вычисление погрешности
      if (eps < maxeps) break;
  }
  ```

• Вычисление и вывод контрольной суммы

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

Для сетки размером [math]\displaystyle{ N \times N }[/math]:

• Каждая итерация требует [math]\displaystyle{ O(N^2) }[/math] операций
• Количество итераций зависит от сходимости, ограничено `itmax`
• Общая сложность: [math]\displaystyle{ O(it\_max \cdot N^2) }[/math]

Конкретно для каждой внутренней точки на каждой итерации:

• 4 операции чтения
• 3 сложения
• 1 деление
• 1 запись
• 1 вычитание и взятие модуля для погрешности

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

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

Информационный граф метода релаксации для [math]\displaystyle{ N }[/math] = 4 , [math]\displaystyle{ max\_it }[/math] = 2

1) Первая группа вершин (входные данные и инициализация) :

Расположена в двумерной области с координатами [math]\displaystyle{ (i, j) }[/math]. Диапазоны изменения координат: [math]\displaystyle{ i ∈ [0, N−1] , j ∈ [0, N−1] }[/math].

Операция:

• если [math]\displaystyle{ i = 0 }[/math] или [math]\displaystyle{ i = N−1 }[/math] или [math]\displaystyle{ j = 0 }[/math] или [math]\displaystyle{ j = N−1 }[/math] → [math]\displaystyle{ A[i][j] = 0 }[/math] (граничные условия),
• иначе [math]\displaystyle{ A[i][j] = 1 + i + j }[/math] (начальное приближение внутри области).

Результат: элементы массива [math]\displaystyle{ A[i][j] }[/math] на итерации 0. Все вершины этой группы являются источниками (входными данными алгоритма). На рисунке соответствуют синим прямоугольникам слева с надписью [math]\displaystyle{ A[i][j]^{(0)} }[/math].

2) Вторая группа вершин (релаксация — вычисление нового слоя) :

Расположена в трёхмерной области с координатами [math]\displaystyle{ (i, j, k) }[/math].

Диапазоны изменения координат: [math]\displaystyle{ i ∈ [1, N−2] , j ∈ [1, N−2] , k ∈ [1, it\_max] }[/math], где [math]\displaystyle{ k }[/math] — номер выполненной итерации.

Соответствующая операция: [math]\displaystyle{ B[i][j] = (A[i−1][j] + A[i+1][j] + A[i][j−1] + A[i][j+1]) / 4.0 }[/math] (усреднение по четырём соседям — шаблон «крест»).

Аргументы операции: четыре значения [math]\displaystyle{ A[⋅][⋅] }[/math] с предыдущего слоя (итерации [math]\displaystyle{ k−1 }[/math]).

Результат срабатывания: элемент нового слоя [math]\displaystyle{ B[i][j]^{(k)} }[/math]. На рисунке — зелёные блоки с надписью [math]\displaystyle{ B[i][j]^{(k)} }[/math].

3) Третья группа вершин (вычисление невязки и обновление [math]\displaystyle{ A }[/math]) :

Расположена в той же трёхмерной области с координатами [math]\displaystyle{ (i, j, k) }[/math]. Операция выполняется сразу после второй группы на той же итерации [math]\displaystyle{ k }[/math] :

1. [math]\displaystyle{ e = |A[i][j]^{(k-1)} − B[i][j]^{(k)}| }[/math]
2. [math]\displaystyle{ A[i][j]^{(k)} = B[i][j]^{(k)} }[/math]
3. [math]\displaystyle{ eps^{(k)} = max(eps^{(k)}, e) }[/math]

Аргументы: значения [math]\displaystyle{ A[i][j]^{(k-1)} }[/math] из первой группы вершин и [math]\displaystyle{ B[i][j]^{(k)} }[/math] из второй группы вершин.

Результат:

• обновлённое значение [math]\displaystyle{ A[i][j]^{(k)} }[/math], которое будет использовано на следующей итерации,
• вклад в глобальную норму невязки [math]\displaystyle{ eps^{(k)} }[/math].

На рисунке — жёлтые блоки с надписями вида [math]\displaystyle{ A[i][j]^{(k)} = B[i][j]^{(k)} }[/math].

4) Четвёртая группа вершин (условия завершения) :

Состоит из двух типов вершин скалярной размерности (координата [math]\displaystyle{ k }[/math]):

• Вершины сравнения [math]\displaystyle{ eps^{(k)} \lt max\_eps }[/math] :
  • Координата [math]\displaystyle{ k ∈ [1, it\_max] }[/math].
  • Аргумент — накопленное на итерации [math]\displaystyle{ k }[/math] значение [math]\displaystyle{ eps^{(k)} }[/math].
  • На рисунке — фиолетовые овалы с надписью [math]\displaystyle{ eps^{(k)} \lt max\_eps }[/math].

• Вершина проверки числа итераций [math]\displaystyle{ iteration = it\_max }[/math] : На рисунке — красный прямоугольник справа.

Обе вершины являются стоками-условиями (управляющие вершины). При выполнении любого из условий (YES) происходит выход из цикла (break). При невыполнении обоих условий (NO) алгоритм переходит к следующей итерации [math]\displaystyle{ k+1 }[/math].

Таким образом, граф имеет ярко выраженную волновую структуру по итерациям [math]\displaystyle{ k }[/math], с полным барьером синхронизации в конце каждой итерации (все элементы нового слоя должны быть вычислены перед подсчётом невязки и переходом к следующей итерации).

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

Алгоритм обладает высоким потенциалом параллелизма:

• На каждой итерации все внутренние точки могут обновляться независимо
• Высота ЯПФ: [math]\displaystyle{ O(it\_max) }[/math] (последовательные итерации)
• Ширина ЯПФ: [math]\displaystyle{ O(N^2) }[/math] (все точки могут обрабатываться параллельно на одной итерации)

В MPI-реализации:

• Сетка разбивается на горизонтальные полосы
• Каждый процесс обрабатывает свою полосу
• Необходим обмен граничными строками между соседними процессами
• Используются неблокирующие операции MPI для перекрытия вычислений и коммуникаций

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

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

• Размер сетки [math]\displaystyle{ N }[/math] (целое число)
• Матрица значений [math]\displaystyle{ A }[/math] размером [math]\displaystyle{ N \times N }[/math] (вещественные числа)
• Параметры: maxeps (точность), itmax (максимальное число итераций)

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

• Матрица [math]\displaystyle{ A }[/math] с решением (вещественные числа)
• Контрольная сумма [math]\displaystyle{ S }[/math] (вещественное число)
• Число выполненных итераций (целое)
• Время выполнения (в параллельной версии)

Объём данных: [math]\displaystyle{ O(N^2) }[/math]

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

1. Детерминированность: При фиксированных входных данных и последовательном выполнении результат детерминирован.
2. Сходимость: Метод сходится медленно (линейно), количество итераций пропорционально [math]\displaystyle{ N^2 }[/math] для достижения заданной точности.
3. Устойчивость: Алгоритм устойчив к ошибкам округления.
4. Масштабируемость: Хорошо масштабируется на параллельных системах благодаря локальности данных.
5. Вычислительная мощность: Отношение числа операций к объёму данных составляет [math]\displaystyle{ O(it\_max) }[/math].

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

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

В последовательной реализации:

1. Используются две матрицы [math]\displaystyle{ A }[/math] и [math]\displaystyle{ B }[/math] для хранения значений на текущей и предыдущей итерациях.
2. Граничные условия заданы нулевыми.
3. Начальное приближение: [math]\displaystyle{ A[i][j] = 1.0 + i + j }[/math] для внутренних точек.
4. Реализован режим вычисления максимальной погрешности для контроля сходимости.
5. Контрольная сумма вычисляется по формуле, включающей веса [math]\displaystyle{ (i+1)(j+1) }[/math].

Ключевые особенности:

   ```c
   // Вычисление нового значения
   B[i][j] = (A[i-1][j] + A[i+1][j] + A[i][j-1] + A[i][j+1]) / 4.0;
   // Обновление и вычисление погрешности
   e = fabs(A[i][j] - B[i][j]);
   A[i][j] = B[i][j];
   eps = Max(eps, e);
   ```

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

MPI-реализация включает следующие особенности:

• Разбиение данных: Сетка разбивается на горизонтальные полосы, каждая обрабатывается отдельным процессом.
• Распределение нагрузки: При неравномерном разбиении (когда [math]\displaystyle{ N }[/math] не кратно числу процессов) первые процессы получают на одну строку больше.
• Обмен граничными данными: Используются неблокирующие операции [math]\displaystyle{ MPI\_Isend / MPI\_Irecv }[/math] для перекрытия вычислений и коммуникаций.
• Вычисление погрешности: Локальные погрешности собираются с помощью [math]\displaystyle{ MPI\_Allreduce }[/math] для нахождения глобального максимума.
• Вычисление контрольной суммы: Частичные суммы вычисляются на каждом процессе и собираются с помощью [math]\displaystyle{ MPI\_Reduce }[/math].

   ```c
   // Обмен граничными строками
   if (rank > 0) {
       MPI_Irecv(localA[0], N, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &req[req_count++]);
       MPI_Isend(localA[1], N, MPI_FLOAT, rank-1, 0, MPI_COMM_WORLD, &req[req_count++]);
   }
   if (rank < size-1) {
       MPI_Irecv(localA[local_rows+1], N, MPI_FLOAT, rank+1, 0, MPI_COMM_WORLD, &req[req_count++]);
       MPI_Isend(localA[local_rows], N, MPI_FLOAT, rank+1, 0, MPI_COMM_WORLD, &req[req_count++]);
   }
   ```

• Оптимизация: Пропуск вычислений для граничных строк на первом и последнем процессах.
• Измерение времени: Использование [math]\displaystyle{ MPI\_Wtime }[/math] для точного измерения времени выполнения.

Результаты прогонов

Для получения результатов была проведена серия экспериментов с различными:

1. Размерами сетки ( [math]\displaystyle{ N }[/math] = 512, 1024, 2048, 4096 )
2. Количеством процессов ( 1, 2, 3, ... , 28 )
3. Параметрами точности ( [math]\displaystyle{ max\_eps }[/math] = 0.1e-7 )

Зависимость времени выполнения от размера матрицы и числа процессов

Результаты экспериментов подтверждают корректность и эффективность MPI-реализации :

• Время выполнения растёт пропорционально квадрату размера сетки, демонстрируя сложность [math]\displaystyle{ O(N²) }[/math].
• Алгоритм хорошо масштабируется: для [math]\displaystyle{ N }[/math] = 4096 использование 28 процессов ускоряет вычисления в 20 раз по сравнению с последовательной версией.
• Значения контрольной суммы сохраняются при любом количестве процессов, что подтверждает корректность обмена данными.
• Наилучшая эффективность достигается при балансе между объёмом вычислений на процесс и затратами на коммуникации.

Литература

1. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. — М.: Наука, 1978.
2. Самарский А. А., Гулин А. В. Численные методы. - M.: Наука, 1989.
3. Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes: The Art of Scientific Computing. — 3rd ed. — Cambridge University Press, 2007.
4. Gropp W., Lusk E., Skjellum A. Using MPI: Portable Parallel Programming with the Message-Passing Interface. — 3rd ed. — MIT Press, 2014.
5. Документация по MPI: https://www.mpi-forum.org/docs/