Перевод FFTW3 этого примера преобразования на месте FFTW2 реальных данных. ⇐ C++

[Гость]

В документации FFTW2 приведен пример прямого и обратного преобразования реальных данных с использованием MPI. Я пытаюсь воссоздать это в FFTW3.

Пример FFTW2:

#include int main(int argc, char **argv) { const int nx = ..., ny = ..., nz = ...; int local_nx, local_x_start, local_ny_after_transpose, local_y_start_after_transpose, total_local_size; интервал х, у, z; план rfftwnd_mpi_plan, iplan; fftw_real *данные, *работа; fftw_complex *cdata; MPI_Init(&argc,&argv); /* создаем планы вперед и назад: */ план = rfftw3d_mpi_create_plan(MPI_COMM_WORLD, Нью-Йорк, Нью-Йорк, Новая Зеландия, ФФТВ_РЕАЛ_ТО_КОМПЛЕКС, ФФТВ_ESTIMATE); iplan = rfftw3d_mpi_create_plan(MPI_COMM_WORLD, /* размеры РЕАЛЬНЫХ данных --> */ nx, ny, nz, ФФТВ_КОМПЛЕКС_ТО_РЕАЛ, ФФТВ_ESTIMATE); rfftwnd_mpi_local_sizes(plan, &local_nx, &local_x_start, &local_ny_after_transpose, &local_y_start_after_transpose, &total_local_size); данные = (fftw_real*) malloc(sizeof(fftw_real) * total_local_size); /* рабочая область имеет тот же размер, что и данные: */ работа = (fftw_real*) malloc(sizeof(fftw_real) * total_local_size); /* инициализируем данные как f(x,y,z): */ for (x = 0; x < local_nx; ++x) для (y = 0; y < ny; ++y) для (z = 0; z < nz; ++z) данные[(x*ny + y) * (2*(nz/2+1)) + z] = f(x + local_x_start, y, z); /* Теперь вычислим прямое преобразование: */ rfftwnd_mpi (план, 1, данные, работа, FFTW_TRANSPOSED_ORDER); /* данные теперь сложные, поэтому приведите указатель к типу: */ cdata = (fftw_complex*) данные; /* умножаем мнимую часть на 2, ради интереса: (обратите внимание, что данные транспонируются) */ for (y = 0; y < local_ny_after_transpose; ++y) для (х = 0; х < nx; ++x) для (z = 0; z < (nz/2+1); ++z) cdata[(y*nx + x) * (nz/2+1) + z].im *= 2,0; /* Наконец, вычисляем обратное преобразование; результат транспонируется обратно в исходное расположение данных: */ rfftwnd_mpi (iplan, 1, данные, работа, FFTW_TRANSPOSED_ORDER); бесплатно (данные); бесплатно (работа); rfftwnd_mpi_destroy_plan (план); rfftwnd_mpi_destroy_plan (iplan); MPI_Finalize(); } Я знаю, что требуется много изменений (и было бы здорово, если бы кто-нибудь где-нибудь перевел эту функцию). Это нормально.

Мои основные вопросы:
[*]Каковы правильные функции замены для rfftw3d_mpi_create_plan и rfftwnd_mpi? Кажется, это могут быть fftw_mpi_plan_dft_r2c_3d и fftw_execute, но у меня возникли трудности, потому что fftw_mpi_plan_dft_r2c_3d требует аргументов double *in и fftw_complex *out. Но поскольку я хочу преобразовать реальные данные на месте, я не могу использовать реальные данные как *in, так и *out, поскольку типы не совпадают. [*]Ни один из примеров с FFTW3, которые я нашел, похоже, не использует переменную рабочей области work. Это больше не требуется?
Я подозреваю, что переведенный код будет выглядеть примерно так

#include int main(int argc, char **argv) { const int nx = ..., ny = ..., nz = ...; ptrdiff_t local_nx, local_x_start, local_ny_after_transpose, local_y_start_after_transpose, total_local_size; интервал х, у, z; план fftw_plan, iplan; fftw_real *данные, *работа; fftw_complex *cdata; MPI_Init(&argc,&argv); /* Находим размеры массива для каждого процессора */ alloc_local = fftw_mpi_local_size_3d_transpose(nx, ny, nz, MPI_COMM_WORLD, &local_nx, &local_x_start, &local_ny_after_transpose, &local_y_start_after_transpose); /* создаем планы вперед и назад: */ план = rfftw3d_mpi_create_plan(MPI_COMM_WORLD, Нью-Йорк, Нью-Йорк, Новая Зеландия, ФФТВ_РЕАЛ_ТО_КОМПЛЕКС, ФФТВ_ESTIMATE); iplan = rfftw3d_mpi_create_plan(MPI_COMM_WORLD, /* размеры РЕАЛЬНЫХ данных --> */ nx, ny, nz, ФФТВ_КОМПЛЕКС_ТО_РЕАЛ, ФФТВ_ESTIMATE); данные = (fftw_real*) malloc(sizeof(fftw_real) * total_local_size); /* Я не знаю, требуется ли/полезна ли рабочая область в FFTW3 */ работа = (fftw_real*) malloc(sizeof(fftw_real) * total_local_size); /* инициализируем данные как f(x,y,z): */ for (x = 0; x < local_nx; ++x) для (y = 0; y < ny; ++y) для (z = 0; z < nz; ++z) данные[(x*ny + y) * (2*(nz/2+1)) + z] = f(x + local_x_start, y, z); /* Код в конечном итоге дает сбой, я думаю, может быть, потому, что я здесь неправильно привел типы? */ plan = fftw_mpi_plan_dft_r2c_3d(nx, ny, nz, data, (fftw_complex*) данные, MPI_COMM_WORLD, FFTW_MPI_TRANSPOSED_OUT); iplan = fftw_mpi_plan_dft_r2c_3d(nx, ny, nz, (fftw_complex*) данные, данные, MPI_COMM_WORLD, FFTW_MPI_TRANSPOSED_IN); /* Теперь вычислим прямое преобразование: */ /* Здесь кодируйте ошибки */ fftw_execute (план); /* данные теперь сложные, поэтому приведите указатель к типу: */ cdata = (fftw_complex*) данные; /* умножаем мнимую часть на 2, ради интереса: (обратите внимание, что данные транспонируются) */ for (y = 0; y < local_ny_after_transpose; ++y) для (х = 0; х < nx; ++x) для (z = 0; z < (nz/2+1); ++z) cdata[(y*nx + x) * (nz/2+1) + z].im *= 2,0; /* Наконец, вычисляем обратное преобразование; результат транспонируется обратно в исходное расположение данных: */ fftw_execute (iplan); бесплатно (данные); бесплатно (работа); fftw_destroy_plan (план); fftw_destroy_plan (iplan); MPI_Finalize(); } Но я понимаю, что это неправильно.