Уравнение Пуассона, решение дискретным преобразованием Фурье

Материал из Алговики
Перейти к навигации Перейти к поиску

Основные авторы описания: В.М.Степаненко, Е.В.Мортиков, Вад.В.Воеводин (раздел 2.2)

Содержание

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

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

Уравнение Пуассона для многомерного пространства имеет вид:

[math] \sum_{i=1}^{N}\frac{\partial^2 \phi}{\partial x_i^2}=f,~\mathbf{x}\in D, [/math]

где [math]D \in \mathbb{R}^N[/math] - область определения решения [math]\phi(\mathbf{x})[/math], [math]\mathbf{x}=(x_1,...,x_N)^T[/math] - вектор независимых переменных. Уравнение Пуассона дополняется граничными условиями:

[math] B(\phi)=F, \mathbf{x} \in \Gamma(D), [/math]

где [math]\Gamma(D)[/math] - граница области [math]D[/math], a [math]B(\phi)[/math] - оператор, определяющий граничные условия. [math]B(\phi)=\phi[/math] соответствует граничным условиям Дирихле, [math]B(\phi)=\partial\phi/\partial n[/math] ([math]\mathbf{n}[/math] - внешняя нормаль к границе [math]\Gamma(D)[/math]) - условиям Неймана. Иногда задают также смешанное граничное условие: [math]B(\phi)=C\phi+\partial\phi/\partial n[/math] ([math]C[/math] - константа). Встречаются также так называемые "периодические граничные условия", при которых задача ставится для бесконечной области, но предполагается периодичность решения по подмножеству переменных из [math]\mathbf{x}[/math].

Уравнение Пуассона возникает во многих задачах математической физики, например, в электростатике (в этом случае [math]\phi[/math] - потенциал электрической силы) и гидродинамике ([math]\phi[/math] - давление жидкости или газа); при этом [math]N=2,3[/math] для плоской и трехмерной задач, соответственно.

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

[math] \sum_{i=1}^{N}\frac{\phi_{k_1,...,k_i+1,...,k_N}-2\phi_{k_1,...,k_i,...,k_N}+\phi_{k_1,...,k_i-1,...,k_N}}{\Delta x_i^2}=f_{k_1,...,k_N},~(k_1,...,k_N) \in D_N. [/math]

Здесь вторые производные заменены конечно-разностной аппроксимацией второго порядка точности (образуя т.н. схему "крест" для плоской задачи), а решение ищется на дискретном множестве точек [math]N[/math]-мерного пространства, [math]D_N[/math]. Конечными разностями аппроксимируются при этом также граничные условия.

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

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

[math] \frac{\phi_{i+1,j,k}-2\phi_{i,j,k}+\phi_{i-1,j,k}}{ \Delta x^2}\,+\, \frac{\phi_{i,j+1,k}-2\phi_{i,j,k}+\phi_{i,j-1,k}}{ \Delta y^2}\,+\, \frac{\phi_{i,j,k+1}-2\phi_{i,j,k}+\phi_{i,j,k-1}}{ \Delta z^2} = f_{i,j,k},~(i,j,k) \in D_h, [/math]

где [math]D_h=\{0:N_x-1\}\times\{0:N_y-1\}\times\{0:N_z-1\} [/math] - параллелепипед в сеточной области. В качестве граничных условий для простоты примем так называемые троякопериодические граничные условия

[math] \begin{align} \phi_{0,j,k} &= \phi_{N_x,j,k},\\ \phi_{i,0,k} &= \phi_{i,N_y,k},\\ \phi_{i,j,0} &= \phi_{i,j,N_z}. \end{align} [/math]

Периодические граничные условия удовлетворяются "автоматически", если представить решение в виде стандартного дискретного обратного преобразования Фурье:

[math] \phi_{i,j,k}=\frac{1}{N_x N_y N_z}\sum_{l=0}^{N_x-1}\sum_{m=0}^{N_y-1}\sum_{n=0}^{N_z-1}\Phi_{l,m,n}e^{2\pi \overline{i} \left(\frac{il}{N_x}+\frac{jm}{N_y}+\frac{kn}{N_z}\right)}. [/math]

Здесь [math]\overline{i}=\sqrt{-1}[/math]. Аналогичное разложение применяется также к правой части уравнения, [math]f_{i,j,k}[/math]. Решение дискретного уравнения Пуассона преобразованием Фурье удобно тем, что базисные функции разложения Фурье являются собственными функциями дискретного оператора Лапласа. А именно, подставляя разложение Фурье неизвестной функции [math]\phi_{i,j,k}[/math] и правой части в исходное уравнение, получаем:

[math] \Phi_{l,m,n}=-\frac{F_{l,m,n}}{4\left[\sin^2\left(\frac{\pi l}{N_x}\right) + \sin^2\left(\frac{\pi m}{N_y}\right) + \sin^2\left(\frac{\pi n}{N_z}\right) \right]}, [/math]

где [math]F_{l,m,n}[/math] - преобразование Фурье правой части. Отсюда очевиден алгоритм решения уравнения: сначала правая часть уравнения разлагается в ряд Фурье, затем по приведенной выше формуле вычисляются коэффициенты Фурье решения и, наконец, решение восстанавливается обратным преобразованием Фурье.

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

Вычислительным ядром алгоритма является одномерное преобразование Фурье. В самом деле, обратное дискретное преобразование Фурье может быть записано как

[math] \phi_{i,j,k}=\frac{1}{N_x} \sum_{l=0}^{N_x-1} \left[ e^{2\pi \overline{i}\left(\frac{il}{N_x}\right)} \frac{1}{N_y} \sum_{m=0}^{N_y-1} \left[ e^{2\pi \overline{i}\left(\frac{jm}{N_y}\right)} \frac{1}{N_z} \sum_{n=0}^{N_z-1} \Phi_{l,m,n}e^{2\pi \overline{i} \left(\frac{kn}{N_z}\right)}\right]\right], [/math]

откуда видно, что трехмерное преобразование Фурье сводится в последовательности трех одномерных. Для вычисления одномерного преобразования широко используется эффективный алгоритм быстрого преобразования Фурье (БПФ, FFT, Fast Fourier Transform)[1].

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

Из приведенного выше ясно, что основной макрооперацией алгоритма решения уравнения Пуассона является одномерное быстрое преобразование Фурье. В дальнейшем будем обозначать ее [math]\text{FFT}_i,~i=x,y,z[/math] по каждому из трех направлений, соответственно, а обратное преобразование Фурье - [math]\text{IFFT}_i,~i=x,y,z[/math].

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

Введем для краткости следующие обозначения сеточных функций: [math]\phi_h=\{\phi_{i,j,k},~(i,j,k) \in D_h\},~F_h=\{F_{l,m,n},~l=0,...,N_x-1,~m=0,...,N_y-1,~n=0,...,N_z-1\}[/math] и аналогично [math]f_h[/math] и [math]\Phi_h[/math]. Тогда алгоритм запишется следующим образом:

  1. Вычисление [math]\text{FFT}_x(f_h)[/math]
  2. Вычисление [math]\text{FFT}_y(\text{FFT}_x(f_h))[/math]
  3. Вычисление [math]\text{FFT}_z(\text{FFT}_y(\text{FFT}_x(f_h)))=F_h[/math]
  4. Расчет [math]\Phi_h[/math] по [math]F_h[/math]
  5. Вычисление [math]\text{IFFT}_x(\Phi_h)[/math]
  6. Вычисление [math]\text{IFFT}_y(\text{IFFT}_x(\Phi_h))[/math]
  7. Вычисление [math]\text{IFFT}_z(\text{IFFT}_y(\text{IFFT}_x(\Phi_h)))=\phi_h[/math]

Примечательно, что весь описанный алгоритм реализуется с помощью одного трехмерного массива [math]N_x \times N_y \times N_z[/math], поскольку результат каждого одномерного преобразования Фурье можно записать во входной массив и затем использовать его на входе для следующего преобразования, а операция (4) - это поэлементное изменение массива.

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

Рассмотрим случай куба, [math]N_x=N_y=N_z=N=2^k[/math]. Тогда сложность одномерного БПФ по каждому из направлений составляет (по порядку) [math]N(\text{log}_2 N)[/math] операций. Поскольку одномерное БПФ на каждом из этапов алгоритма, кроме четвертого, осуществляется [math]N^2[/math] раз, а количество операций на этапе 4 составляет (снова по порядку) [math]N^3[/math], то общее количество арифметических действий составляет [math]N^3(6\text{log}_2 N+1)[/math].

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

Представим информационный граф алгоритма для этапов 1-4, поскольку этапы 5-7 реализуются аналогично этапам 1-3.

Зеленые прямоугольники обозначают одномерные сечения входного трехмерного массива. Красным кругом обозначена операция одномерного БПФ, например [math]FFT_{x11}[/math] означает БПФ вдоль оси x при [math]j=1,~k=1[/math]. Желтыми кругами представлены поэлементные преобразования массива для получения коэффициентов Фурье решения (этап 4 алгоритма), например [math]L_{2N,1,...,N}[/math] означает преобразование элементов [math]i=2,~j=N,~k=1,...,N[/math]. Для удобства целочисленные индексы нумеруются числами [math]1,...,N[/math] вместо [math]0,...,N-1[/math]

Информационная зависимость между ярусами графа заключается в том, что одномерное БПФ зависит от результата произведенного ранее БПФ по перпендикулярному направлению (т.е. с предыдущего яруса), только если обрабатываемые этими БПФ одномерные сечения массива пересекаются. Следовательно, каждое из БПФ по второй координате, [math]FFT_{i,y,k},~i=1,...,N[/math], зависит от результатов каждого из БПФ по первой координате [math]FFT_{x,j,k},~j=1,...,N[/math], и только от них, если эти преобразования производятся в одной плоскости [math]k=K[/math]. Соответственно, каждое БПФ по третьей координате [math]FFT_{i,j,z},~j=1,...,N[/math], зависит от каждого БПФ по второй координате [math]FFT_{i,y,k},~k=1,...,N[/math], в одной плоскости [math]i=I[/math].

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

В описанном выше алгоритме можно выделить, по меньшей мере, два уровня параллелизма. Во-первых, вычисление каждого одномерного БПФ может быть распределено между вычислительными ядрами. И, во-вторых, на каждом этапе алгоритма (ярусе графа) одномерные БПФ независимы друг от друга и также могут быть выполнены паралельно (координатный параллелизм). Для сопоставления параллельной (при координатном распараллеливании) и последовательной сложности алгоритма определим в качестве основных операций одномерное БПФ и поэлементный перерасчет массива (этап 4 алгоритма). Тогда последовательный алгоритм будет состоять из последовательно выполняемых [math]6N^2[/math] одномерных БПФ и [math]N^3[/math] операций над отдельными элементами массива, параллельный же алгоритм будет выполнен за 6 шагов одномерного БПФ и 1 операцию поэлементного преобразования. Если же выразить введенные макрооперации через элементарные арифметические операции, то получим, что сложность последовательного алгоритма - [math]6N^3(\text{log}_2N)+N^3[/math], а сложность параллельного - [math]6N(\text{log}_2N)+1[/math]. Тогда при больших [math]N[/math] отношение последовательной сложности к параллельной становится равным [math]\propto N^2(6\text{log}_2N \text{ln}2+\text{ln}2+5)[/math].

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

В качестве входных данных выступает трехмерный массив [math]N \times N \times N[/math] правой части уравнения. Выходной массив - трехмерный массив решения той же размерности.

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

Как было показано выше, при неограниченном количестве вычислительных устройств и больших объемах входных данных вычислительная сложность параллельного алгоритма убывает быстрее, чем [math]N^2[/math].

Поскольку объем входных данных составляет [math]N^3[/math], вычислительная мощность последовательного алгоритма становится равной [math]6\text{log}_2N+1[/math], а параллельного алгоритма - [math](6N\text{log}_2N+1)/N^3 \propto N^{-2}[/math] (при больших [math]N[/math]).

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

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

Приведем пример реализации алгоритма на языке Фортран-90 с применением библиотеки FFTW. В этом случае трехмерное прямое преобразование Фурье правой части уравнения реализуется как последовательность трех одномерных (этапы 1-3 в п.1.5). Полученные коэффициенты Фурье делятся на собственное число оператора Лапласа (этап 4 в п.1.5) (см. комментарии в коде): 

         ! x - transform
         do k = 1, k1
           do j = 1, j1
             workfor_x(1:i1) = cmplx(f_r(1:i1,j,k),f_i(1:i1,j,k))
             call FFTW_EXECUTE_DFT(planfor_x,workfor_x,workfor_x)
             f_r(1:i1,j,k) = real(real (workfor_x(1:i1)) )
             f_i(1:i1,j,k) = real(aimag(workfor_x(1:i1)) )
           enddo
         enddo
         ! y - transform
         do k = 1, k1
           do i = 1, i1
             workfor_y(1:j1) = cmplx(f_r(i,1:j1,k),f_i(i,1:j1,k))
             call FFTW_EXECUTE_DFT(planfor_y,workfor_y,workfor_y)
             f_r(i,1:j1,k) = real(real (workfor_y(1:j1)) )
             f_i(i,1:j1,k) = real(aimag(workfor_y(1:j1)) )
           enddo
         enddo
         ! z - transform and dividing by eigenvalue
         do j = 1, j1
           do i = 1, i1
             workfor_z(1:k1) = cmplx(f_r(i,j,1:k1),f_i(i,j,1:k1))
             call FFTW_EXECUTE_DFT(planfor_z,workfor_z,workfor_z)
             f_r(i,j,1:k1) = real(real (workfor_z(1:k1)) )
             f_i(i,j,1:k1) = real(aimag(workfor_z(1:k1)) )
             ! Dividing by eigenvalue of Laplace operator
             sinx = sin((i-1) * Pi / i1)^2 * 4
             siny = sin((j-1) * Pi / j1)^2 * 4
             do k = 1, k1
               sinz = sin((k-1) * Pi / k1)^2 * 4
               f_r(i,j,1:k1) = f_r(i,j,1:k1)/(sinx + siny + sinz)
               f_i(i,j,1:k1) = f_i(i,j,1:k1)/(sinx + siny + sinz)
             enddo
           enddo
         enddo

Как только коэффициенты Фурье решения найдены, аналогично производится обратное преобразование Фурье: 

         ! inverse x - transform
         do k = 1, k1
           do j = 1, j1
             workback_x(1:i1) = cmplx(f_r(1:i1,j,k),f_i(1:i1,j,k))
             call FFTW_EXECUTE_DFT(planback_x,workback_x,workback_x)
             f_r(1:i1,j,k) = real(real (workback_x(1:i1)) )/real(i1)
             f_i(1:i1,j,k) = real(aimag(workback_x(1:i1)) )/real(i1)
           enddo
         enddo
         ! inverse y - transform
         do k = 1, k1
           do i = 1, i1
             workback_y(1:j1) = cmplx(f_r(i,1:j1,k),f_i(i,1:j1,k))
             call FFTW_EXECUTE_DFT(planback_y,workback_y,workback_y)
             f_r(i,1:j1,k) = real(real (workback_y(1:j1)) )/real(j1)
             f_i(i,1:j1,k) = real(aimag(workback_y(1:j1)) )/real(j1)
           enddo
         enddo
         ! inverse z - transform
         do j = 1, j1
           do i = 1, i1
             workback_z(1:k1) = cmplx(f_r(i,j,1:k1),f_i(i,j,1:k1))
             call FFTW_EXECUTE_DFT(planback_z,workback_z,workback_z)
             f_r(i,j,1:k1) = real(real (workback_z(1:k1)) )
             f_i(i,j,1:k1) = real(aimag(workback_z(1:k1)) )
           enddo
         enddo

Как видим, вложенность циклов соответствует порядку расположения элементов массивов в языке Фортран - по столбцам - это позволяет существенно сократить время обращения к памяти. Многократно используются одни и те же вспомогательные массивы, что при небольших размерах сетки должно обеспечить быстроту работы с кэш-памятью. Однако, данная реализация не является оптимальной. Так, вместо многократного вызова подпрограммы FFTW_EXECUTE_DFT в циклах, можно было бы обращаться к одной подпрограмме, выполняющей заданное множество одномерных БПФ.

Из приведенных выше листингов программы становится очевиден также и ресурс параллелизма алгоритма: каждое одномерное преобразование Фурье выполняется независимо. Самым простым способом использования этого обстоятельства будет вставка соответствующих директив OpenMP перед внутренними циклами - в случае использования многоядерной архтитектуры компьютера на общей памяти. В случае же параллельного выполнения на кластерах с распределенной памятью, необходимо будет позаботиться о распределении входного массива (правой части уравнения) между процессорами, а также "транспозиции" массивов при переходе от БПФ по одному из направлений к БПФ по следующему направлению.

Заметим также, что в приведенной выше реализации для расчета одномерного БПФ производится вызов процедуры библиотеки FFTW, обрабатывающей массив комплексных чисел. В то же время, существует вариант БПФ для преобразования действительных массивов, ускоряющий БПФ в два раза (более подробно см. на сайте проекта FFTW).

Важным обстоятельством является то, что большинство вариантов БПФ, включая алгоритм Кули-Тьюки, в арифметике с плавающей запятой имеют очень низкую погрешность. Верхняя граница относительной ошибки для алгоритма Кули-Тьюки составляет [math]O(\epsilon \log N)[/math] (для сравнения - относительная ошибка преобразования Фурье непосредственно по формуле дискретного преобразования Фурье - [math]O(\epsilon N^{3/2})[/math]) [2]. Поэтому для ускорения алгоритма можно использовать одинарную точность представления переменных без существенного увеличения ошибки.

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

2.2.1 Локальность реализации алгоритма

2.2.1.1 Структура обращений в память и качественная оценка локальности
Рисунок 3. Уравнение Пуассона, решение дискретным преобразованием Фурье. Общий профиль обращений в память

На рис. 3 представлен профиль обращений в память для реализации уравнения Пуассона, решение дискретным преобразованием Фурье. Данный профиль состоит из обращений к нескольким служебным массивам (фрагмент 1) и основному массиву (франмент 2). Видно, что профиль делится на два практически одинаковых этапа (вертикальная желтая линия). В целом можно отметить, что общее количество элементов невелико (чуть более 1000 элементов), при этом профиль составляет около 100 тысяч обращений. При этом видно, что в верхней части профиля плотность обращений выше, что говорит, вероятно, о высокой локальности вычислений; однако обращения к основному массиву достаточно разрознены. Перейдем к более детальному рассмотрению локальности.

На рис. 4 представлен фрагмент 1, на котором расположены все обращения к служебным массивам. Исходя только из данного графика, нельзя сказать, какова структура обращения внутри блоков, однако в данном случае это и не требуется. Общее число задействованных адресов мало – всего порядка 100 элементов, при этом они используются в течение всей программы. Более того, они организованы в блочную структуру, что улучшает работу с памятью. Поэтому можно с уверенностью сказать, что локальность (как пространственная, так и временная) в данном случае очень высока.

Рисунок 4. Профиль обращений, фрагмент 1

Перейдем к рассмотрению основного массива. На рис. 5 представлен фрагмент 2 от общего профиля, показанного на рис. 3. Можно сразу же отметить следующий факт. На графике представлено максимум несколько сотен обращений, а по оси Х стоит отметка более 35000. Это значит, что значительно больше обращений происходит вне этого фрагмента, а именно в служебных массивах. То есть в целом обращения к основному массиву происходят не так часто. Это подтверждается при рассмотрении исходного кода программы: данные из основного массива на самом деле только копируются в служебные массивы, где над ними выполняются преобразование Фурье, после чего они копируются обратно.

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

Рисунок 5. Профиль обращений, фрагмент 2

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

2.2.1.2 Количественная оценка локальности

Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен здесь (функция Kernel). Условия запуска описаны здесь.

Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.

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

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

Рисунок 6. Сравнение значений оценки daps

Вторая характеристика – cvg – предназначена для получения более машинно-независимой оценки локальности. Она определяет, насколько часто в программе необходимо подгружать данные в кэш-память. Соответственно, чем меньше значение cvg, тем реже это нужно делать, тем лучше локальность.

На рис. 7 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем, в общем случае, выше локальность). Можно увидеть, что, как и в случае алгоритма Кули-Тьюки, локальность оценивается как очень высокая. Поскольку в данной программе реализации этого алгоритма посвящена большая часть общего профиля, то и причины такого расхождения между локальностью и производительностью в данных случаях одинаковы (см. подробное рассмотрение алгоритма Кули-Тьюки).

Рисунок 7. Сравнение значений оценки cvg

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

Как было показано выше, существует, по меньшей мере, два уровня параллелизма в описываемом алгоритме. Во-первых, это параллельная реализация одномерных БПФ и, во-вторых, координатный параллелизм - параллельное исполнение независимых одномерных БПФ [3]. Эти два уровня могут быть скомбинированы, например, путем распределения разных БПФ между MPI-процессами, а каждого отдельного БПФ - между нитями на общей памяти (POSIX, OpenMP, и т.д.). Ниже мы приведем детали базового алгоритма распараллеливания с использованием координатного параллелизма, рассмотрим его недостатки и перечислим возможные способы повышения его масштабируемости.

Координатный параллелизм в контексте данной задачи реализуется в результате разбиения трехмерной области на подобласти, в каждой из которых можно независимо от расчетов в других подобластях, выполнить одномерные БПФ. К настоящему времени предложены две реализации такого подхода: одномерное разбиение области (slab decomposition) и двумерное разбиение (pencil decomposition) [4]. При одномерном разбиении требуется меньшее количество MPI-пересылок, однако, максимальное количество процессоров/ядер ограничено максимальным по трем измерениям размером области. Поэтому с точки зрения массивного параллелизма двумерная декомпозиция предпочтительнее, и ниже мы остановимся на этом подходе.

При MPI-реализации необходимо принять во внимание накладные расходы, вызванные пересылками данных между процессорами. Для реализации обсуждаемого алгоритма необходимо осуществлять т.н. транспонирования массивов, при которых на каждом из трех этапов преобразования Фурье (этапы 1-3 и 5-7 в п.1.5) обеспечивается нахождение элементов массива вдоль соответствующих направлений в памяти соответствующих процессоров:

Транспонирование трехмерного массива при параллельной реализации трехмерного БПФ

Рассмотрим соотношение вычислений и пересылок данных при MPI-реализации. Вычислительная сложность алгоритма при покоординатном распараллеливании, при запуске на [math]M[/math] MPI-процессах составляет [math] \left(6N^3(\log_2 N) + N^3\right)/M[/math], т.е., в лучшем случае, [math] 6N(\log_2 N) + N[/math] (на [math]N^2[/math] процессах). При этом, суммарное количество отсылаемых и получаемых MPI-пакетов составляет [math]8N_{px}+4N_{py}[/math], а суммарное количество отсылаемых и получаемых элементов массива на каждом процессе - [math]12N^3/(N_{px}N_{py})[/math] (здесь - [math]N_{px},~N_{py}[/math] - количества процессов в двумерной MPI-декомпозиции области по осям [math]x[/math] и [math]y[/math], соответственно - см. рис. выше). Учитывая, что [math]M=N_{px}N_{py}[/math], получаем, что отношение количества пересылаемых элементов массива к количеству арифметических операций составляет [math]\alpha=12/(6\log_2 N + 1)[/math]. Таким образом, отношение объема пересылаемых элементов к количеству операций не зависит от количества MPI-процессов и уменьшается с ростом размера задачи. При больших [math]N[/math], однако, эта зависимость становится слабой, т.е. "насыщается". В то же время, при фиксированном размере задачи, рост количества используемых процессов, не изменяя параметра [math]\alpha[/math], вызывает рост числа пакетов ([math]\propto \sqrt{M}[/math], в силу [math]N_{px} \propto N_{py} \propto \sqrt{M}[/math] и формул, приведенных выше). В силу латентности сети, это приводит к росту времени пересылки единицы данных и, следовательно, замедлению ускорения алгоритма.

Приведенный здесь анализ может служить только для ориентировочных оценок, поскольку он не дает времен вычислений и MPI-обменов. Более детальные расчеты для случаев одномерной и двумерной декомпозиции области обширно представлены в литературе [5][6].

В профилируемом в следующих разделах программном коде параллельной реализации алгоритма вызовы одномерных БПФ по трем направлениям, показанные в листинге п.2.1, перемежаются вызовами специальных процедур, осуществляющих транспонирование массива. Отдельно транспонируются действительная и мнимая части массива (что уже не оптимально - лучше пересылать комплексные числа) с помощью коммуникаций "точка-точка". При этом используются асинхронные пересылки с помощью MPI_ISEND и MPI_IRECV.

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

Набор и границы значений изменяемых параметров запуска реализации алгоритма:

  • число процессоров [4 : 100] со значениями квадрата целого числа;
  • размер области [100x100x100 : 250x100x100] с шагом 50x100x100.

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

  • минимальная эффективность реализации 0,132%;
  • максимальная эффективность реализации 1,34%.

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

Рисунок 8. Параллельная реализация алгоритма. Изменение производительности в зависимости от числа процессоров и размера области.
Рисунок 9. Параллельная реализация алгоритма. Изменение эффективности в зависимости от числа процессоров и размера области.

Приведенные результаты хорошо качественно согласуются с оценками, приведенными в предыдущем разделе. Так, на рис.8 видно, что с ростом [math]N[/math] масштабируемость алгоритма повышается - производительность алгоритма быстрее растет с ростом числа ядер. Это согласуется с уменьшением [math]\alpha[/math] при увеличении [math]N[/math], хотя выражение [math]\alpha=12/(6\log_2 N+1)[/math] дает значительно более слабую зависимость. В целом же, алгоритм масштабируется достаточно слабо: при максимальном из задаваемых размеров области 250х100х100 при росте числа ядер в 25 раз (от 4 до 100), производительность вычислений вырастает только в 6 раз.

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

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

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

Ниже представлены известные к настоящему моменту реализации трехмерного БПФ и их основные характеристики.

Реализация Разбиение области Технологии распараллеливания Обмены в MPI
FFTW 1D MPI+OpenMP
FFTE 1D & 2D MPI+OpenMP MPI_ALLTOALL
P3DFFT 1D & 2D MPI
cuFFT
AccFFT 1D & 2D MPI, CUDA MPI_ALLTOALL
MKL FFT 1D MPI
PFFT MPI

3 Литература

  1. ↑ Марчук Г.И. Методы вычислительной математики, М., 1977, 456с.
  2. ↑ Gentleman, W. M.; Sande, G. (1966). "Fast Fourier transforms—for fun and profit". Proc. AFIPS 29: 563–578. doi:10.1145/1464291.1464352
  3. ↑ Anshu Dubey and Daniele Tessera. Redistribution strategies for portable parallel FFT: a case study. Concurrency and Computation: Practice and Experience, 13(3):209–220, 2001.
  4. ↑ Orlando Ayala, Lian-Ping Wang, Parallel implementation and scalability analysis of 3D Fast Fourier Transform using 2D domain decomposition, Parallel Computing, Volume 39, Issue 1, January 2013, Pages 58-77, ISSN 0167-8191, http://dx.doi.org/10.1016/j.parco.2012.12.002.
  5. ↑ Orlando Ayala, Lian-Ping Wang, Parallel implementation and scalability analysis of 3D Fast Fourier Transform using 2D domain decomposition, Parallel Computing, Volume 39, Issue 1, January 2013, Pages 58-77, ISSN 0167-8191, http://dx.doi.org/10.1016/j.parco.2012.12.002.
  6. ↑ P. Dmitruk, L.-P. Wang, W.H. Matthaeus, R. Zhang, D. Seckel, Scalable parallel FFT for spectral simulations on a Beowulf cluster, Parallel Comput. 27 (2001) 1921–1936.

Основные авторы описания: В.М.Степаненко, Е.В.Мортиков

Источник — https://algowiki-project.org/w/ru/index.php?title=Уравнение_Пуассона,_решение_дискретным_преобразованием_Фурье&oldid=3365