Chris@10: Chris@10:
Chris@10:Chris@10: Next: Wisdom of Fortran?, Chris@10: Previous: FFTW Execution in Fortran, Chris@10: Up: Calling FFTW from Legacy Fortran Chris@10:
In C, you might have something like the following to transform a Chris@10: one-dimensional complex array: Chris@10: Chris@10:
fftw_complex in[N], out[N]; Chris@10: fftw_plan plan; Chris@10: Chris@10: plan = fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_ESTIMATE); Chris@10: fftw_execute(plan); Chris@10: fftw_destroy_plan(plan); Chris@10:Chris@10:
In Fortran, you would use the following to accomplish the same thing: Chris@10: Chris@10:
double complex in, out Chris@10: dimension in(N), out(N) Chris@10: integer*8 plan Chris@10: Chris@10: call dfftw_plan_dft_1d(plan,N,in,out,FFTW_FORWARD,FFTW_ESTIMATE) Chris@10: call dfftw_execute_dft(plan, in, out) Chris@10: call dfftw_destroy_plan(plan) Chris@10:Chris@10:
Chris@10: Notice how all routines are called as Fortran subroutines, and the
Chris@10: plan is returned via the first argument to dfftw_plan_dft_1d
.
Chris@10: Notice also that we changed fftw_execute
to
Chris@10: dfftw_execute_dft
(see FFTW Execution in Fortran). To do
Chris@10: the same thing, but using 8 threads in parallel (see Multi-threaded FFTW), you would simply prefix these calls with:
Chris@10:
Chris@10:
integer iret Chris@10: call dfftw_init_threads(iret) Chris@10: call dfftw_plan_with_nthreads(8) Chris@10:Chris@10:
Chris@10: (You might want to check the value of iret
: if it is zero, it
Chris@10: indicates an unlikely error during thread initialization.)
Chris@10:
Chris@10:
To transform a three-dimensional array in-place with C, you might do: Chris@10: Chris@10:
fftw_complex arr[L][M][N]; Chris@10: fftw_plan plan; Chris@10: Chris@10: plan = fftw_plan_dft_3d(L,M,N, arr,arr, Chris@10: FFTW_FORWARD, FFTW_ESTIMATE); Chris@10: fftw_execute(plan); Chris@10: fftw_destroy_plan(plan); Chris@10:Chris@10:
In Fortran, you would use this instead: Chris@10: Chris@10:
double complex arr Chris@10: dimension arr(L,M,N) Chris@10: integer*8 plan Chris@10: Chris@10: call dfftw_plan_dft_3d(plan, L,M,N, arr,arr, Chris@10: & FFTW_FORWARD, FFTW_ESTIMATE) Chris@10: call dfftw_execute_dft(plan, arr, arr) Chris@10: call dfftw_destroy_plan(plan) Chris@10:Chris@10:
Chris@10: Note that we pass the array dimensions in the “natural” order in both C Chris@10: and Fortran. Chris@10: Chris@10:
To transform a one-dimensional real array in Fortran, you might do: Chris@10: Chris@10:
double precision in Chris@10: dimension in(N) Chris@10: double complex out Chris@10: dimension out(N/2 + 1) Chris@10: integer*8 plan Chris@10: Chris@10: call dfftw_plan_dft_r2c_1d(plan,N,in,out,FFTW_ESTIMATE) Chris@10: call dfftw_execute_dft_r2c(plan, in, out) Chris@10: call dfftw_destroy_plan(plan) Chris@10:Chris@10:
Chris@10: To transform a two-dimensional real array, out of place, you might use Chris@10: the following: Chris@10: Chris@10:
double precision in Chris@10: dimension in(M,N) Chris@10: double complex out Chris@10: dimension out(M/2 + 1, N) Chris@10: integer*8 plan Chris@10: Chris@10: call dfftw_plan_dft_r2c_2d(plan,M,N,in,out,FFTW_ESTIMATE) Chris@10: call dfftw_execute_dft_r2c(plan, in, out) Chris@10: call dfftw_destroy_plan(plan) Chris@10:Chris@10:
Chris@10: Important: Notice that it is the first dimension of the Chris@10: complex output array that is cut in half in Fortran, rather than the Chris@10: last dimension as in C. This is a consequence of the interface routines Chris@10: reversing the order of the array dimensions passed to FFTW so that the Chris@10: Fortran program can use its ordinary column-major order. Chris@10: Chris@10: Chris@10: Chris@10: Chris@10: