PARALLEL.RU - Информационно-аналитический центр по параллельным вычислениям

КЛАСТЕР

Библиотека FFTW для выполнения быстрых преобразований Фурье (БПФ)

• Сайт разработчиков библиотеки - www.fftw.org.
  • Документация к версии 2.1.4 библиотеки.
• Смотрите также: сервер по численному анализу НИВЦ МГУ.

---

Введение

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

1. Одномерное преобразование Фурье для комплексных чисел
2. Многомерное преобразование Фурье для комплексных чисел
3. Одномерное преобразование Фурье для действительных чисел
4. Многомерное преобразование Фурье для действительных чисел

Работа с FFTW

Сначала происходит построение "плана", который оптимизирует время вычислений для данной задачи. Затем построенный план передается в качестве параметра функциям, которые непосредственно отвечают за вычисление БПФ.

Одномерное преобразование Фурье для комплексных чисел:

#include<fftw.h>
#define N 100

int main(void)
{
int i;
fftw_complex in[N], out[N];
fftw_plan plan,plan_inv;

/* создать оценочный план для прямого БПФ на массиве из N точек */
plan = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);

/* создать оценочный план для обратного БПФ на массиве из N точек */
plan_inv = fftw_create_plan(N, FFTW_BACKWARD, FFTW_ESTIMATE);

/* вычислить спектр сигнала, заданного массивом in, в соответствии с планом plan */
fftw_one(plan, in, out);

/* вычислительный код программы

...

*/
/* вычислить сигнала по спектру, заданному массивом out, в соответствии с планом plan_inv */
fftw_one(plan_inv, out, in);

/* после обратного преобразования Фурье, полученные данные следует нормировать на длину исходного массива */
for(i=0; i < N; i++)
{
in[i] /= N;
}
/* удаление ранее созданного плана */
fftw_destroy_plan(plan);
fftw_destroy_plan(plan_inv);
return 0;
}

Пояснения к функции fftw_create_plan

fftw_create_plan(int size, fftw_direction dir, int flags)

• size - количество точек, используемых для вычисления БПФ.
• dir - флаг прямого FFTW_FORWARD или обратного FFTW_BACKWARD преобразования Фурье.
• flag - FFTW_ESTIMATE создает план, который не является полностью оптимальным для данной конфигурации; FFTW_MEASURE - используется для создания оптимального плана, в этом случае значительно увеличивается время работы по вычислению плана и БПФ.

Также можно использовать флаг FFTW_IN_PLACE:

	plan = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE)

тогда в качестве out указывается значение NULL:

	fftw_one(plan, in, NULL)

Компиляция под Unix

Как правило, для сборки исполняемого модуля нужно указать а) директорию, где находятся заголовочные файлы (параметр -I); б) директорию, где находятся библиотечные файлы (параметр -L); в) один или несколько библиотечных файлов (параметр -l); Например:

	mpicc -c -I/usr/local/fftw/include offtw.c 
	mpicc offtw.o -L/usr/local/fftw/lib -lfftw -lm

Многомерные преобразования для комплексных чисел

#include<fftw.h>
#define N 100
#define M 100
#define K 100

int main(void)
{
int i;
fftw_complex *in;
fftwnd_plan plan,plan_inv;

in=(fftw_complex*)malloc((K*(N*M+M)+K)*sizeof(fftw_complex));

/* создать оценочный план для прямого БПФ на массиве из N точек */
plan = fftw3d_create_plan(N, M, K, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);

/* создать оценочный план для обратного БПФ на массиве из N точек */
plan_inv = fftw3d_create_plan(N, M, K, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);

/* вычислить спектр сигнала, заданного массивом in, в соответствии с планом plan */
fftwnd_one(plan, in, NULL);

/* вычислительный код программы

...

*/
/* вычислить сигнала по спектру, заданному массивом in, в соответствии с планом plan_inv */
fftwnd_one(plan_inv, in, NULL);

/* после обратного преобразования Фурье, полученные данные следует нормировать на длину исходной области */
for(i=0; i < K*(N*M+M)+K; i++)
{
	in[i].re/= N*M*K;
	in[i].im/= N*M*K;

}

/* удаление ранее созданного плана */
fftwnd_destroy_plan(plan);
fftwnd_destroy_plan(plan_inv);
return 0;
}

Для построения многомерного плана (вместо приставки fftw_ используется fftwnd_, где n размерность) используются функции:

• fftwnd_plan fftw2d_create_plan(int nx, int ny, fftw_direction dir, int flags) - для двумерного БПФ;
• fftwnd_plan fftw3d_create_plan(int nx, int ny, int nz, fftw_direction dir, int flags) - для трехмерного БПФ.

Здесь nx, ny, nz - количество точек по соответствующим декартовым координатам (x, y, z); Флаги dir и flag имеют те же значения, что и в случае одномерного БПФ.

Функция преобразования сигнал/спектр:

void fftwnd_one(fftwnd_plan plan, fftw_complex *in, fftw_complex *out);

Функция удаления созданного плана:

fftwnd_destroy_plan(fftwnd_plan plan)

Преобразование Фурье для действительных чисел

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

	mpicc -I/usr/local/fftw/include rfftw_test.c -L/usr/local/fftw/lib -lrfftw -lfftw -lm

Модуль параллельной обработки БПФ для MPI

Этот модуль позволяет работать на многопроцессорных машинах с общей и распределенной памятью. Основное его отличие от однопроцессорного модуля FFTW заключается в том, что на каждом процессоре обрабатывается свое подмножество точек. Функция создания плана для параллельной обработки принимает, в дополнение к стандартным параметрам, значение коммуникатора для набора процессов (MPI_COMM_WORLD или другой коммуникатор). Обмен данными производиться с помощью функций MPI_Alltoall или MPI_Alltoallv в зависимости от алгоритма распределения данных по процессорам. Надо отметить, что наилучшее распределение по процессорам происходит при условии, что число точек кратно P², где P - число процессов.

В исходный текст программы нужно включить заголовочный файл в дополнение к . При сборке программы нужно добавить специальную библиотеку с помощью флага -lfftw_mpi, в дополнение к общей библиотеке FFTW (флаг -lfftw).

Пример выполнения одномерного комплексного БПФ с помощью параллельной версии FFTW:

#include<stdlib.h>
#include<mpi.h>
#include<fftw_mpi.h>
#include<fftw.h>
#include<math.h>

int main(int argc,char **argv)
{
int size,k,i,rank,Nx=1000;
int local_n,local_start,local_n_after_transform,local_start_after_transform,total_local_size;
double time;
fftw_complex *A;
fftw_mpi_plan plan,pinv;

MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&size);
MPI_Comm_rank(MPI_COMM_WORLD,&rank);

/* создать план для прямого БПФ, в MPI версии FFTW указывать флаг FFTW_IN_PLACE не нужно*/
plan = fftw_mpi_create_plan(MPI_COMM_WORLD, Nx, FFTW_FORWARD, FFTW_ESTIMATE);
/* функция определяющая длину и порядок данных на конкретном процессором */
fftw_mpi_local_sizes(plan, &local_n,&local_start, &local_n_after_transform, &local_start_after_transform, &total_local_size);

A=(fftw_complex*)malloc(total_local_size*sizeof(fftw_complex));
/* обнуление массива
:.
*/
/* заполнение массива данными для дальнейшего счета */
for(k=local_start,i=0;k <local_start+local_n;k++,i++)
{
A[i].re = f(k);
A[i].im = f(k);
}

/* выполнение прямого БПФ */
fftw_mpi(plan,1,A,NULL);


pinv = fftw_mpi_create_plan(MPI_COMM_WORLD, Nx, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_mpi_local_sizes(pinv, &local_n, &local_start, &local_n_after_transform, &local_start_aftеr_transform, &total_local_size);

fftw_mpi(pinv,1,A,NULL);
for(i=0;i < total_local_size;i++)
{
A[i].re/=Nx;
A[i].im/=Nx;
}
fftw_mpi_destroy_plan(plan).
fftw_mpi_destroy_plan(pinv).

MPI_Finalize();
return 0;
}

Комментарии к используемым функциям

Функция создания плана для одномерного комплексного БПФ:

	fftw_mpi_plan fftw_mpi_create_plan(MPI_Comm comm, int nx, fftw_direction dir, int flags);

Параметры:

• comm - коммуникатор группы процессоров, на которых будет производиться вычисление БПФ;
• nx - длина массива данных на всех процессорах;
• dir - флаг прямого FFTW_FORWARD или обратного FFTW_BACKWARD преобразования Фурье.
• flag - FFTW_ESTIMATE создает план, который не является полностью оптимальным для данной конфигурации; FFTW_MEASURE - используется для создания оптимального плана, в этом случае значительно увеличивается время работы по вычислению плана и БПФ. FFTW_IN_PLACE не указывается, т.к. в этом случае он стоит по умолчанию. Если не важен порядок выходных данных можно использовать FFTW_SCRAMBLER_OUTPUT для прямого и FFTW_SCRAMBLER_INPUT для обратного БПФ.

Вычисление количества точек, обрабатываемых на каждом из процессоров: void fftw_mpi_local_sizes(fftw_mpi_plan plan,int *local_n,int *local_start, int *local_n_after_transform, int *local_start_after_transform, int *total_local_size);

• plan - ранее созданный план
• local_n - количество точек на процессоре перед преобразованием;
• local_start - указывает с какой точки начинается нумерация на данном процессоре перед преобразованием; (таким образом, на каждом процессоре перед преобразованием размещаются точки от local_start до local_start+local_n-1)
• local_n_after_transform - количество точек на процессоре после преобразования;
• local_start_after_transform - указывает с какой точки начинается нумерация на данном процессоре после преобразования;
• total_local_size - размер рабочей области на процессоре

Функция параллельного выполнения БПФ:

 
	fftw_mpi(fftw_mpi_plan pinv, n_fields, fftw_complex *local_data, fftw_complex *work)

Здесь n_fields - количество одинаковых векторов, которые надо преобразовать (в простейшем случае равен 1). Для выполнения многомерных преобразований в параллельной версии используются аналогичные функции с префиксом fftwnd_mpi.

---

© Лаборатория Параллельных Информационных Технологий, НИВЦ МГУ