Distributed-memory FFTW with MPI
********************************

In this chapter we document the parallel FFTW routines for parallel systems supporting the MPI message-passing interface.  Unlike the shared-memory threads described in the previous chapter, MPI allows you to use distributed-memory parallelism, where each CPU has its own separate memory, and which can scale up to clusters of many thousands of processors.  This capability comes at a price, however: each process only stores a portion of the data to be transformed, which means that the data structures and programming-interface are quite different from the serial or threads versions of FFTW.

Distributed-memory parallelism is especially useful when you are transforming arrays so large that they do not fit into the memory of a single processor.  The storage per-process required by FFTW's MPI routines is proportional to the total array size divided by the number of processes.  Conversely, distributed-memory parallelism can easily pose an unacceptably high communications overhead for small problems; the threshold problem size for which parallelism becomes advantageous will depend on the precise problem you are interested in, your hardware, and your MPI implementation.

A note on terminology: in MPI, you divide the data among a set of ``processes'' which each run in their own memory address space. Generally, each process runs on a different physical processor, but this is not required.  A set of processes in MPI is described by an opaque data structure called a ``communicator,'' the most common of which is the predefined communicator MPI_COMM_WORLD which refers to all processes.  For more information on these and other concepts common to all MPI programs, we refer the reader to the documentation at the MPI home page (http://www.mcs.anl.gov/research/projects/mpi/).

We assume in this chapter that the reader is familiar with the usage of the serial (uniprocessor) FFTW, and focus only on the concepts new to the MPI interface.

FFTW MPI Installation
=====================

All of the FFTW MPI code is located in the mpi subdirectory of the FFTW package.  On Unix systems, the FFTW MPI libraries and header files are automatically configured, compiled, and installed along with the uniprocessor FFTW libraries simply by including --enable-mpi in the flags to the configure script (see Installation on Unix).

Any implementation of the MPI standard, version 1 or later, should work with FFTW.  The configure script will attempt to automatically detect how to compile and link code using your MPI implementation.  In some cases, especially if you have multiple different MPI implementations installed or have an unusual MPI software package, you may need to provide this information explicitly.

Most commonly, one compiles MPI code by invoking a special compiler command, typically mpicc for C code.  The configure script knows the most common names for this command, but you can specify the MPI compilation command explicitly by setting the MPICC variable, as in ./configure MPICC=mpicc ....

If, instead of a special compiler command, you need to link a certain library, you can specify the link command via the MPILIBS variable, as in ./configure MPILIBS=-lmpi ....  Note that if your MPI library is installed in a non-standard location (one the compiler does not know about by default), you may also have to specify the location of the library and header files via LDFLAGS and CPPFLAGS variables, respectively, as in ./configure LDFLAGS=-L/path/to/mpi/libs CPPFLAGS=-I/path/to/mpi/include ....

Linking and Initializing MPI FFTW
=================================

Programs using the MPI FFTW routines should be linked with -lfftw3_mpi -lfftw3 -lm on Unix in double precision, -lfftw3f_mpi -lfftw3f -lm in single precision, and so on (see Precision). You will also need to link with whatever library is responsible for MPI on your system; in most MPI implementations, there is a special compiler alias named mpicc to compile and link MPI code.

Before calling any FFTW routines except possibly fftw_init_threads (see Combining MPI and Threads), but after calling MPI_Init, you should call the function:

     void fftw_mpi_init(void);

If, at the end of your program, you want to get rid of all memory and other resources allocated internally by FFTW, for both the serial and MPI routines, you can call:

     void fftw_mpi_cleanup(void);

which is much like the fftw_cleanup() function except that it also gets rid of FFTW's MPI-related data.  You must not execute any previously created plans after calling this function.

2d MPI example
==============

Before we document the FFTW MPI interface in detail, we begin with a simple example outlining how one would perform a two-dimensional N0 by N1 complex DFT.

     #include <fftw3-mpi.h>

     int main(int argc, char **argv)
     {
         const ptrdiff_t N0 = ..., N1 = ...;
         fftw_plan plan;
         fftw_complex *data;
         ptrdiff_t alloc_local, local_n0, local_0_start, i, j;

         MPI_Init(&argc, &argv);
         fftw_mpi_init();

         /* get local data size and allocate */
         alloc_local = fftw_mpi_local_size_2d(N0, N1, MPI_COMM_WORLD,
                                              &local_n0, &local_0_start);
         data = fftw_alloc_complex(alloc_local);

         /* create plan for in-place forward DFT */
         plan = fftw_mpi_plan_dft_2d(N0, N1, data, data, MPI_COMM_WORLD,
                                     FFTW_FORWARD, FFTW_ESTIMATE);

         /* initialize data to some function my_function(x,y) */
         for (i = 0; i < local_n0; ++i) for (j = 0; j < N1; ++j)
            data[i*N1 + j] = my_function(local_0_start + i, j);

         /* compute transforms, in-place, as many times as desired */
         fftw_execute(plan);

         fftw_destroy_plan(plan);

         MPI_Finalize();
     }

As can be seen above, the MPI interface follows the same basic style of allocate/plan/execute/destroy as the serial FFTW routines.  All of the MPI-specific routines are prefixed with fftw_mpi_ instead of fftw_.  There are a few important differences, however:

First, we must call fftw_mpi_init() after calling MPI_Init (required in all MPI programs) and before calling any other fftw_mpi_ routine.

Second, when we create the plan with fftw_mpi_plan_dft_2d, analogous to fftw_plan_dft_2d, we pass an additional argument: the communicator, indicating which processes will participate in the transform (here MPI_COMM_WORLD, indicating all processes). Whenever you create, execute, or destroy a plan for an MPI transform, you must call the corresponding FFTW routine on all processes in the communicator for that transform.  (That is, these are collective calls.)  Note that the plan for the MPI transform uses the standard fftw_execute and fftw_destroy routines (on the other hand, there are MPI-specific new-array execute functions documented below).

Third, all of the FFTW MPI routines take ptrdiff_t arguments instead of int as for the serial FFTW.  ptrdiff_t is a standard C integer type which is (at least) 32 bits wide on a 32-bit machine and 64 bits wide on a 64-bit machine.  This is to make it easy to specify very large parallel transforms on a 64-bit machine.  (You can specify 64-bit transform sizes in the serial FFTW, too, but only by using the guru64 planner interface.  See 64-bit Guru Interface.)

Fourth, and most importantly, you don't allocate the entire two-dimensional array on each process.  Instead, you call fftw_mpi_local_size_2d to find out what portion of the array resides on each processor, and how much space to allocate. Here, the portion of the array on each process is a local_n0 by N1 slice of the total array, starting at index local_0_start.  The total number of fftw_complex numbers to allocate is given by the alloc_local return value, which may be greater than local_n0 * N1 (in case some intermediate calculations require additional storage).  The data distribution in FFTW's MPI interface is described in more detail by the next section.

Given the portion of the array that resides on the local process, it is straightforward to initialize the data (here to a function myfunction) and otherwise manipulate it.  Of course, at the end of the program you may want to output the data somehow, but synchronizing this output is up to you and is beyond the scope of this manual.  (One good way to output a large multi-dimensional distributed array in MPI to a portable binary file is to use the free HDF5 library; see the HDF home page (http://www.hdfgroup.org/).)

MPI Data Distribution
=====================

The most important concept to understand in using FFTW's MPI interface is the data distribution.  With a serial or multithreaded FFT, all of the inputs and outputs are stored as a single contiguous chunk of memory.  With a distributed-memory FFT, the inputs and outputs are broken into disjoint blocks, one per process.

In particular, FFTW uses a 1d block distribution of the data, distributed along the first dimension.  For example, if you want to perform a 100,200 complex DFT, distributed over 4 processes, each process will get a 25,200 slice of the data. That is, process 0 will get rows 0 through 24, process 1 will get rows 25 through 49, process 2 will get rows 50 through 74, and process 3 will get rows 75 through 99.  If you take the same array but distribute it over 3 processes, then it is not evenly divisible so the different processes will have unequal chunks.  FFTW's default choice in this case is to assign 34 rows to processes 0 and 1, and 32 rows to process 2.

FFTW provides several fftw_mpi_local_size routines that you can call to find out what portion of an array is stored on the current process.  In most cases, you should use the default block sizes picked by FFTW, but it is also possible to specify your own block size.  For example, with a 100,200 array on three processes, you can tell FFTW to use a block size of 40, which would assign 40 rows to processes 0 and 1, and 20 rows to process 2.  FFTW's default is to divide the data equally among the processes if possible, and as best it can otherwise.  The rows are always assigned in ``rank order,'' i.e. process 0 gets the first block of rows, then process 1, and so on.  (You can change this by using MPI_Comm_split to create a new communicator with re-ordered processes.)  However, you should always call the fftw_mpi_local_size routines, if possible, rather than trying to predict FFTW's distribution choices.

In particular, it is critical that you allocate the storage size that is returned by fftw_mpi_local_size, which is not necessarily the size of the local slice of the array.  The reason is that intermediate steps of FFTW's algorithms involve transposing the array and redistributing the data, so at these intermediate steps FFTW may require more local storage space (albeit always proportional to the total size divided by the number of processes).  The fftw_mpi_local_size functions know how much storage is required for these intermediate steps and tell you the correct amount to allocate.

Basic and advanced distribution interfaces
------------------------------------------

As with the planner interface, the fftw_mpi_local_size distribution interface is broken into basic and advanced (_many) interfaces, where the latter allows you to specify the block size manually and also to request block sizes when computing multiple transforms simultaneously.  These functions are documented more exhaustively by the FFTW MPI Reference, but we summarize the basic ideas here using a couple of two-dimensional examples.

For the 100,200 complex-DFT example, above, we would find the distribution by calling the following function in the basic interface:

     ptrdiff_t fftw_mpi_local_size_2d(ptrdiff_t n0, ptrdiff_t n1, MPI_Comm comm,
                                      ptrdiff_t *local_n0, ptrdiff_t *local_0_start);

Given the total size of the data to be transformed (here, n0 = 100 and n1 = 200) and an MPI communicator (comm), this function provides three numbers.

First, it describes the shape of the local data: the current process should store a local_n0 by n1 slice of the overall dataset, in row-major order (n1 dimension contiguous), starting at index local_0_start.  That is, if the total dataset is viewed as a n0 by n1 matrix, the current process should store the rows local_0_start to local_0_start+local_n0-1.  Obviously, if you are running with only a single MPI process, that process will store the entire array: local_0_start will be zero and local_n0 will be n0.  See Row-major Format.

Second, the return value is the total number of data elements (e.g., complex numbers for a complex DFT) that should be allocated for the input and output arrays on the current process (ideally with fftw_malloc or an fftw_alloc function, to ensure optimal alignment).  It might seem that this should always be equal to local_n0 * n1, but this is not the case.  FFTW's distributed FFT algorithms require data redistributions at intermediate stages of the transform, and in some circumstances this may require slightly larger local storage.  This is discussed in more detail below, under Load balancing.

The advanced-interface local_size function for multidimensional transforms returns the same three things (local_n0, local_0_start, and the total number of elements to allocate), but takes more inputs:

     ptrdiff_t fftw_mpi_local_size_many(int rnk, const ptrdiff_t *n,
                                        ptrdiff_t howmany,
                                        ptrdiff_t block0,
                                        MPI_Comm comm,
                                        ptrdiff_t *local_n0,
                                        ptrdiff_t *local_0_start);

The two-dimensional case above corresponds to rnk = 2 and an array n of length 2 with n[0] = n0 and n[1] = n1. This routine is for any rnk > 1; one-dimensional transforms have their own interface because they work slightly differently, as discussed below.

First, the advanced interface allows you to perform multiple transforms at once, of interleaved data, as specified by the howmany parameter.  (hoamany is 1 for a single transform.)

Second, here you can specify your desired block size in the n0 dimension, block0.  To use FFTW's default block size, pass FFTW_MPI_DEFAULT_BLOCK (0) for block0.  Otherwise, on P processes, FFTW will return local_n0 equal to block0 on the first P / block0 processes (rounded down), return local_n0 equal to n0 - block0 * (P / block0) on the next process, and local_n0 equal to zero on any remaining processes.  In general, we recommend using the default block size (which corresponds to n0 / P, rounded up).

For example, suppose you have P = 4 processes and n0 = 21.  The default will be a block size of 6, which will give local_n0 = 6 on the first three processes and local_n0 = 3 on the last process.  Instead, however, you could specify block0 = 5 if you wanted, which would give local_n0 = 5 on processes 0 to 2, local_n0 = 6 on process 3.  (This choice, while it may look superficially more ``balanced,'' has the same critical path as FFTW's default but requires more communications.)

Load balancing
--------------

Ideally, when you parallelize a transform over some P processes, each process should end up with work that takes equal time. Otherwise, all of the processes end up waiting on whichever process is slowest.  This goal is known as ``load balancing.''  In this section, we describe the circumstances under which FFTW is able to load-balance well, and in particular how you should choose your transform size in order to load balance.

Load balancing is especially difficult when you are parallelizing over heterogeneous machines; for example, if one of your processors is a old 486 and another is a Pentium IV, obviously you should give the Pentium more work to do than the 486 since the latter is much slower. FFTW does not deal with this problem, however---it assumes that your processes run on hardware of comparable speed, and that the goal is therefore to divide the problem as equally as possible.

For a multi-dimensional complex DFT, FFTW can divide the problem equally among the processes if: (i) the first dimension n0 is divisible by P; and (ii), the product of the subsequent dimensions is divisible by P.  (For the advanced interface, where you can specify multiple simultaneous transforms via some ``vector'' length howmany, a factor of howmany is included in the product of the subsequent dimensions.)

For a one-dimensional complex DFT, the length N of the data should be divisible by P squared to be able to divide the problem equally among the processes.

Transposed distributions
------------------------

Internally, FFTW's MPI transform algorithms work by first computing transforms of the data local to each process, then by globally transposing the data in some fashion to redistribute the data among the processes, transforming the new data local to each process, and transposing back.  For example, a two-dimensional n0 by n1 array, distributed across the n0 dimension, is transformd by: (i) transforming the n1 dimension, which are local to each process; (ii) transposing to an n1 by n0 array, distributed across the n1 dimension; (iii) transforming the n0 dimension, which is now local to each process; (iv) transposing back.

However, in many applications it is acceptable to compute a multidimensional DFT whose results are produced in transposed order (e.g., n1 by n0 in two dimensions).  This provides a significant performance advantage, because it means that the final transposition step can be omitted.  FFTW supports this optimization, which you specify by passing the flag FFTW_MPI_TRANSPOSED_OUT to the planner routines.  To compute the inverse transform of transposed output, you specify FFTW_MPI_TRANSPOSED_IN to tell it that the input is transposed.  In this section, we explain how to interpret the output format of such a transform.

Suppose you have are transforming multi-dimensional data with (at least two) dimensions ndims.  As always, it is distributed along the first dimension 0.  Now, if we compute its DFT with the FFTW_MPI_TRANSPOSED_OUT flag, the resulting output data are stored with the first two dimensions transposed: ndimstrans, distributed along the 1 dimension.  Conversely, if we take the ndimstrans data and transform it with the FFTW_MPI_TRANSPOSED_IN flag, then the format goes back to the original ndims array.

There are two ways to find the portion of the transposed array that resides on the current process.  First, you can simply call the appropriate local_size function, passing ndimstrans (the transposed dimensions).  This would mean calling the local_size function twice, once for the transposed and once for the non-transposed dimensions.  Alternatively, you can call one of the local_size_transposed functions, which returns both the non-transposed and transposed data distribution from a single call. For example, for a 3d transform with transposed output (or input), you might call:

     ptrdiff_t fftw_mpi_local_size_3d_transposed(
                     ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t n2, MPI_Comm comm,
                     ptrdiff_t *local_n0, ptrdiff_t *local_0_start,
                     ptrdiff_t *local_n1, ptrdiff_t *local_1_start);

Here, local_n0 and local_0_start give the size and starting index of the n0 dimension for the non-transposed data, as in the previous sections.  For transposed data (e.g. the output for FFTW_MPI_TRANSPOSED_OUT), local_n1 and local_1_start give the size and starting index of the n1 dimension, which is the first dimension of the transposed data (n1 by n0 by n2).

(Note that FFTW_MPI_TRANSPOSED_IN is completely equivalent to performing FFTW_MPI_TRANSPOSED_OUT and passing the first two dimensions to the planner in reverse order, or vice versa.  If you pass both the FFTW_MPI_TRANSPOSED_IN and FFTW_MPI_TRANSPOSED_OUT flags, it is equivalent to swapping the first two dimensions passed to the planner and passing neither flag.)

One-dimensional distributions
-----------------------------

For one-dimensional distributed DFTs using FFTW, matters are slightly more complicated because the data distribution is more closely tied to how the algorithm works.  In particular, you can no longer pass an arbitrary block size and must accept FFTW's default; also, the block sizes may be different for input and output.  Also, the data distribution depends on the flags and transform direction, in order for forward and backward transforms to work correctly.

     ptrdiff_t fftw_mpi_local_size_1d(ptrdiff_t n0, MPI_Comm comm,
                     int sign, unsigned flags,
                     ptrdiff_t *local_ni, ptrdiff_t *local_i_start,
                     ptrdiff_t *local_no, ptrdiff_t *local_o_start);

This function computes the data distribution for a 1d transform of size n0 with the given transform sign and flags. Both input and output data use block distributions.  The input on the current process will consist of local_ni numbers starting at index local_i_start; e.g. if only a single process is used, then local_ni will be n0 and local_i_start will be 0.  Similarly for the output, with local_no numbers starting at index local_o_start.  The return value of fftw_mpi_local_size_1d will be the total number of elements to allocate on the current process (which might be slightly larger than the local size due to intermediate steps in the algorithm).

As mentioned above (see Load balancing), the data will be divided equally among the processes if n0 is divisible by the square of the number of processes.  In this case, local_ni will equal local_no.  Otherwise, they may be different.

For some applications, such as convolutions, the order of the output data is irrelevant.  In this case, performance can be improved by specifying that the output data be stored in an FFTW-defined ``scrambled'' format.  (In particular, this is the analogue of transposed output in the multidimensional case: scrambled output saves a communications step.)  If you pass FFTW_MPI_SCRAMBLED_OUT in the flags, then the output is stored in this (undocumented) scrambled order.  Conversely, to perform the inverse transform of data in scrambled order, pass the FFTW_MPI_SCRAMBLED_IN flag.

In MPI FFTW, only composite sizes n0 can be parallelized; we have not yet implemented a parallel algorithm for large prime sizes.

Multi-dimensional MPI DFTs of Real Data
=======================================

FFTW's MPI interface also supports multi-dimensional DFTs of real data, similar to the serial r2c and c2r interfaces.  (Parallel one-dimensional real-data DFTs are not currently supported; you must use a complex transform and set the imaginary parts of the inputs to zero.)

The key points to understand for r2c and c2r MPI transforms (compared to the MPI complex DFTs or the serial r2c/c2r transforms), are:

  * Just as for serial transforms, r2c/c2r DFTs transform ndims real data to/from ndimshalf complex data: the last dimension of the complex data is cut in half (rounded down), plus one.  As for the serial transforms, the sizes you pass to the plan_dft_r2c and plan_dft_c2r are the ndims dimensions of the real data.

  * Although the real data is conceptually ndims, it is physically stored as an ndimspad array, where the last dimension has been padded to make it the same size as the complex output.  This is much like the in-place serial r2c/c2r interface (see Multi-Dimensional DFTs of Real Data), except that in MPI the padding is required even for out-of-place data.  The extra padding numbers are ignored by FFTW (they are not like zero-padding the transform to a larger size); they are only used to determine the data layout.

  * The data distribution in MPI for both the real and complex data is determined by the shape of the complex data.  That is, you call the appropriate local size function for the ndimshalf complex data, and then use the same distribution for the real data except that the last complex dimension is replaced by a (padded) real dimension of twice the length.

For example suppose we are performing an out-of-place r2c transform of L,M,N real data [padded to L,M,2(N/2+1)], resulting in L,M,N/2+1 complex data.  Similar to the example in 2d MPI example, we might do something like:

     #include <fftw3-mpi.h>

     int main(int argc, char **argv)
     {
         const ptrdiff_t L = ..., M = ..., N = ...;
         fftw_plan plan;
         double *rin;
         fftw_complex *cout;
         ptrdiff_t alloc_local, local_n0, local_0_start, i, j, k;

         MPI_Init(&argc, &argv);
         fftw_mpi_init();

         /* get local data size and allocate */
         alloc_local = fftw_mpi_local_size_3d(L, M, N/2+1, MPI_COMM_WORLD,
                                              &local_n0, &local_0_start);
         rin = fftw_alloc_real(2 * alloc_local);
         cout = fftw_alloc_complex(alloc_local);

         /* create plan for out-of-place r2c DFT */
         plan = fftw_mpi_plan_dft_r2c_3d(L, M, N, rin, cout, MPI_COMM_WORLD,
                                         FFTW_MEASURE);

         /* initialize rin to some function my_func(x,y,z) */
         for (i = 0; i < local_n0; ++i)
            for (j = 0; j < M; ++j)
              for (k = 0; k < N; ++k)
            rin[(i*M + j) * (2*(N/2+1)) + k] = my_func(local_0_start+i, j, k);

         /* compute transforms as many times as desired */
         fftw_execute(plan);

         fftw_destroy_plan(plan);

         MPI_Finalize();
     }

Note that we allocated rin using fftw_alloc_real with an argument of 2 * alloc_local: since alloc_local is the number of complex values to allocate, the number of real values is twice as many.  The rin array is then local_n0,M,2(N/2+1) in row-major order, so its (i,j,k) element is at the index (i*M + j) * (2*(N/2+1)) + k (see Multi-dimensional Array Format).

As for the complex transforms, improved performance can be obtained by specifying that the output is the transpose of the input or vice versa (see Transposed distributions).  In our L,M,N r2c example, including FFTW_TRANSPOSED_OUT in the flags means that the input would be a padded L,M,2(N/2+1) real array distributed over the L dimension, while the output would be a M,L,N/2+1 complex array distributed over the M dimension.  To perform the inverse c2r transform with the same data distributions, you would use the FFTW_TRANSPOSED_IN flag.

Other multi-dimensional Real-Data MPI Transforms
================================================

FFTW's MPI interface also supports multi-dimensional r2r transforms of all kinds supported by the serial interface (e.g. discrete cosine and sine transforms, discrete Hartley transforms, etc.).  Only multi-dimensional r2r transforms, not one-dimensional transforms, are currently parallelized.

These are used much like the multidimensional complex DFTs discussed above, except that the data is real rather than complex, and one needs to pass an r2r transform kind (fftw_r2r_kind) for each dimension as in the serial FFTW (see More DFTs of Real Data).

For example, one might perform a two-dimensional L,M that is an REDFT10 (DCT-II) in the first dimension and an RODFT10 (DST-II) in the second dimension with code like:

         const ptrdiff_t L = ..., M = ...;
         fftw_plan plan;
         double *data;
         ptrdiff_t alloc_local, local_n0, local_0_start, i, j;

         /* get local data size and allocate */
         alloc_local = fftw_mpi_local_size_2d(L, M, MPI_COMM_WORLD,
                                              &local_n0, &local_0_start);
         data = fftw_alloc_real(alloc_local);

         /* create plan for in-place REDFT10 x RODFT10 */
         plan = fftw_mpi_plan_r2r_2d(L, M, data, data, MPI_COMM_WORLD,
                                     FFTW_REDFT10, FFTW_RODFT10, FFTW_MEASURE);

         /* initialize data to some function my_function(x,y) */
         for (i = 0; i < local_n0; ++i) for (j = 0; j < M; ++j)
            data[i*M + j] = my_function(local_0_start + i, j);

         /* compute transforms, in-place, as many times as desired */
         fftw_execute(plan);

         fftw_destroy_plan(plan);

Notice that we use the same local_size functions as we did for complex data, only now we interpret the sizes in terms of real rather than complex values, and correspondingly use fftw_alloc_real.

FFTW MPI Transposes
===================

The FFTW's MPI Fourier transforms rely on one or more global transposition step for their communications.  For example, the multidimensional transforms work by transforming along some dimensions, then transposing to make the first dimension local and transforming that, then transposing back.  Because global transposition of a block-distributed matrix has many other potential uses besides FFTs, FFTW's transpose routines can be called directly, as documented in this section.

Basic distributed-transpose interface
-------------------------------------

In particular, suppose that we have an n0 by n1 array in row-major order, block-distributed across the n0 dimension.  To transpose this into an n1 by n0 array block-distributed across the n1 dimension, we would create a plan by calling the following function:

     fftw_plan fftw_mpi_plan_transpose(ptrdiff_t n0, ptrdiff_t n1,
                                       double *in, double *out,
                                       MPI_Comm comm, unsigned flags);

The input and output arrays (in and out) can be the same.  The transpose is actually executed by calling fftw_execute on the plan, as usual.

The flags are the usual FFTW planner flags, but support two additional flags: FFTW_MPI_TRANSPOSED_OUT and/or FFTW_MPI_TRANSPOSED_IN.  What these flags indicate, for transpose plans, is that the output and/or input, respectively, are locally transposed.  That is, on each process input data is normally stored as a local_n0 by n1 array in row-major order, but for an FFTW_MPI_TRANSPOSED_IN plan the input data is stored as n1 by local_n0 in row-major order.  Similarly, FFTW_MPI_TRANSPOSED_OUT means that the output is n0 by local_n1 instead of local_n1 by n0.

To determine the local size of the array on each process before and after the transpose, as well as the amount of storage that must be allocated, one should call fftw_mpi_local_size_2d_transposed, just as for a 2d DFT as described in the previous section:

     ptrdiff_t fftw_mpi_local_size_2d_transposed
                     (ptrdiff_t n0, ptrdiff_t n1, MPI_Comm comm,
                      ptrdiff_t *local_n0, ptrdiff_t *local_0_start,
                      ptrdiff_t *local_n1, ptrdiff_t *local_1_start);

Again, the return value is the local storage to allocate, which in this case is the number of real (double) values rather than complex numbers as in the previous examples.

Advanced distributed-transpose interface
----------------------------------------

The above routines are for a transpose of a matrix of numbers (of type double), using FFTW's default block sizes.  More generally, one can perform transposes of tuples of numbers, with user-specified block sizes for the input and output:

     fftw_plan fftw_mpi_plan_many_transpose
                     (ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t howmany,
                      ptrdiff_t block0, ptrdiff_t block1,
                      double *in, double *out, MPI_Comm comm, unsigned flags);

In this case, one is transposing an n0 by n1 matrix of howmany-tuples (e.g. howmany = 2 for complex numbers). The input is distributed along the n0 dimension with block size block0, and the n1 by n0 output is distributed along the n1 dimension with block size block1.  If FFTW_MPI_DEFAULT_BLOCK (0) is passed for a block size then FFTW uses its default block size.  To get the local size of the data on each process, you should then call fftw_mpi_local_size_many_transposed.

An improved replacement for MPI_Alltoall
----------------------------------------

We close this section by noting that FFTW's MPI transpose routines can be thought of as a generalization for the MPI_Alltoall function (albeit only for floating-point types), and in some circumstances can function as an improved replacement.

MPI_Alltoall is defined by the MPI standard as:

     int MPI_Alltoall(void *sendbuf, int sendcount, MPI_Datatype sendtype,
                      void *recvbuf, int recvcnt, MPI_Datatype recvtype,
                      MPI_Comm comm);

In particular, for double* arrays in and out, consider the call:

     MPI_Alltoall(in, howmany, MPI_DOUBLE, out, howmany MPI_DOUBLE, comm);

This is completely equivalent to:

     MPI_Comm_size(comm, &P);
     plan = fftw_mpi_plan_many_transpose(P, P, howmany, 1, 1, in, out, comm, FFTW_ESTIMATE);
     fftw_execute(plan);
     fftw_destroy_plan(plan);

That is, computing a P,P transpose on P processes, with a block size of 1, is just a standard all-to-all communication.

However, using the FFTW routine instead of MPI_Alltoall may have certain advantages.  First of all, FFTW's routine can operate in-place (in == out) whereas MPI_Alltoall can only operate out-of-place.

Second, even for out-of-place plans, FFTW's routine may be faster, especially if you need to perform the all-to-all communication many times and can afford to use FFTW_MEASURE or FFTW_PATIENT.  It should certainly be no slower, not including the time to create the plan, since one of the possible algorithms that FFTW uses for an out-of-place transpose is simply to call MPI_Alltoall.  However, FFTW also considers several other possible algorithms that, depending on your MPI implementation and your hardware, may be faster.

FFTW MPI Wisdom
===============

FFTW's ``wisdom'' facility (see Words of Wisdom-Saving Plans) can be used to save MPI plans as well as to save uniprocessor plans. However, for MPI there are several unavoidable complications.

First, the MPI standard does not guarantee that every process can perform file I/O (at least, not using C stdio routines)---in general, we may only assume that process 0 is capable of I/O. (In fact, even this assumption is not technically guaranteed by the standard, although it seems to be universal in actual MPI implementations and is widely assumed by MPI-using software.  Technically, you need to query the MPI_IO attribute of MPI_COMM_WORLD with MPI_Attr_get.  If this attribute is MPI_PROC_NULL, no I/O is possible.  If it is MPI_ANY_SOURCE, any process can perform I/O.  Otherwise, it is the rank of a process that can perform I/O ... but since it is not guaranteed to yield the same rank on all processes, you have to do an MPI_Allreduce of some kind if you want all processes to agree about which is going to do I/O. And even then, the standard only guarantees that this process can perform output, but not input. See e.g. Parallel Programming with MPI by P. S. Pacheco, section 8.1.3.  Needless to say, in our experience virtually no MPI programmers worry about this.) So, if we want to export the wisdom from a single process to a file, we must first export the wisdom to a string, then send it to process 0, then write it to a file.

Second, in principle we may want to have separate wisdom for every process, since in general the processes may run on different hardware even for a single MPI program.  However, in practice FFTW's MPI code is designed for the case of homogeneous hardware (see Load balancing), and in this case it is convenient to use the same wisdom for every process.  Thus, we need a mechanism to synchronize the wisdom.

To address both of these problems, FFTW provides the following two functions:

     void fftw_mpi_broadcast_wisdom(MPI_Comm comm);
     void fftw_mpi_gather_wisdom(MPI_Comm comm);

Given a communicator comm, fftw_mpi_broadcast_wisdom will broadcast the wisdom from process 0 to all other processes. Conversely, fftw_mpi_gather_wisdom will collect wisdom from all processes onto process 0.  (If the plans created for the same problem by different processes are not the same, fftw_mpi_gather_wisdom will arbitrarily choose one of the plans.)  Both of these functions may result in suboptimal plans for different processes if the processes are running on non-identical hardware.  Both of these functions are collective calls, which means that they must be executed by all processes in the communicator.

So, for example, a typical code snippet to import wisdom from a file and use it on all processes would be:

     {
         int rank;

         fftw_mpi_init();
         MPI_Comm_rank(MPI_COMM_WORLD, &rank);
         if (rank == 0) fftw_import_wisdom_from_filename("mywisdom");
         fftw_mpi_broadcast_wisdom(MPI_COMM_WORLD);
     }

(Note that we must call fftw_mpi_init before importing any wisdom that might contain MPI plans.)  Similarly, a typical code snippet to export wisdom from all processes to a file is:

     {
         int rank;

         fftw_mpi_gather_wisdom(MPI_COMM_WORLD);
         MPI_Comm_rank(MPI_COMM_WORLD, &rank);
         if (rank == 0) fftw_export_wisdom_to_filename("mywisdom");
     }

Avoiding MPI Deadlocks
======================

An MPI program can deadlock if one process is waiting for a message from another process that never gets sent.  To avoid deadlocks when using FFTW's MPI routines, it is important to know which functions are collective: that is, which functions must always be called in the same order from every process in a given communicator.  (For example, MPI_Barrier is the canonical example of a collective function in the MPI standard.)

The functions in FFTW that are always collective are: every function beginning with fftw_mpi_plan, as well as fftw_mpi_broadcast_wisdom and fftw_mpi_gather_wisdom. Also, the following functions from the ordinary FFTW interface are collective when they are applied to a plan created by an fftw_mpi_plan function: fftw_execute, fftw_destroy_plan, and fftw_flops.

FFTW MPI Performance Tips
=========================

In this section, we collect a few tips on getting the best performance out of FFTW's MPI transforms.

First, because of the 1d block distribution, FFTW's parallelization is currently limited by the size of the first dimension. (Multidimensional block distributions may be supported by a future version.) More generally, you should ideally arrange the dimensions so that FFTW can divide them equally among the processes. See Load balancing.

Second, if it is not too inconvenient, you should consider working with transposed output for multidimensional plans, as this saves a considerable amount of communications.  See Transposed distributions.

Third, the fastest choices are generally either an in-place transform or an out-of-place transform with the FFTW_DESTROY_INPUT flag (which allows the input array to be used as scratch space).  In-place is especially beneficial if the amount of data per process is large.

Fourth, if you have multiple arrays to transform at once, rather than calling FFTW's MPI transforms several times it usually seems to be faster to interleave the data and use the advanced interface.  (This groups the communications together instead of requiring separate messages for each transform.)

Combining MPI and Threads
=========================

In certain cases, it may be advantageous to combine MPI (distributed-memory) and threads (shared-memory) parallelization. FFTW supports this, with certain caveats.  For example, if you have a cluster of 4-processor shared-memory nodes, you may want to use threads within the nodes and MPI between the nodes, instead of MPI for all parallelization.

In particular, it is possible to seamlessly combine the MPI FFTW routines with the multi-threaded FFTW routines (see Multi-threaded FFTW). However, some care must be taken in the initialization code, which should look something like this:

     int threads_ok;

     int main(int argc, char **argv)
     {
         int provided;
         MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &provided);
         threads_ok = provided >= MPI_THREAD_FUNNELED;

         if (threads_ok) threads_ok = fftw_init_threads();
         fftw_mpi_init();

         ...
         if (threads_ok) fftw_plan_with_nthreads(...);
         ...

         MPI_Finalize();
     }

First, note that instead of calling MPI_Init, you should call MPI_Init_threads, which is the initialization routine defined by the MPI-2 standard to indicate to MPI that your program will be multithreaded.  We pass MPI_THREAD_FUNNELED, which indicates that we will only call MPI routines from the main thread.  (FFTW will launch additional threads internally, but the extra threads will not call MPI code.)  (You may also pass MPI_THREAD_SERIALIZED or MPI_THREAD_MULTIPLE, which requests additional multithreading support from the MPI implementation, but this is not required by FFTW.)  The provided parameter returns what level of threads support is actually supported by your MPI implementation; this must be at least MPI_THREAD_FUNNELED if you want to call the FFTW threads routines, so we define a global variable threads_ok to record this.  You should only call fftw_init_threads or fftw_plan_with_nthreads if threads_ok is true.  For more information on thread safety in MPI, see the MPI and Threads (http://www.mpi-forum.org/docs/mpi-20-html/node162.htm) section of the MPI-2 standard.

Second, we must call fftw_init_threads before fftw_mpi_init.  This is critical for technical reasons having to do with how FFTW initializes its list of algorithms.

Then, if you call fftw_plan_with_nthreads(N), every MPI process will launch (up to) N threads to parallelize its transforms.

For example, in the hypothetical cluster of 4-processor nodes, you might wish to launch only a single MPI process per node, and then call fftw_plan_with_nthreads(4) on each process to use all processors in the nodes.

This may or may not be faster than simply using as many MPI processes as you have processors, however.  On the one hand, using threads within a node eliminates the need for explicit message passing within the node.  On the other hand, FFTW's transpose routines are not multi-threaded, and this means that the communications that do take place will not benefit from parallelization within the node. Moreover, many MPI implementations already have optimizations to exploit shared memory when it is available, so adding the multithreaded FFTW on top of this may be superfluous.

FFTW MPI Reference
==================

This chapter provides a complete reference to all FFTW MPI functions, datatypes, and constants.  See also FFTW Reference for information on functions and types in common with the serial interface.

MPI Files and Data Types
------------------------

All programs using FFTW's MPI support should include its header file:

     #include <fftw3-mpi.h>

Note that this header file includes the serial-FFTW fftw3.h header file, and also the mpi.h header file for MPI, so you need not include those files separately.

You must also link to both the FFTW MPI library and to the serial FFTW library.  On Unix, this means adding -lfftw3_mpi -lfftw3 -lm at the end of the link command.

Different precisions are handled as in the serial interface: see Precision.  That is, fftw_ functions become fftwf_ (in single precision) etcetera, and the libraries become -lfftw3f_mpi -lfftw3f -lm etcetera on Unix.  Long-double precision is supported in MPI, but quad precision (fftwq_) is not due to the lack of MPI support for this type.

MPI Initialization
------------------

Before calling any other FFTW MPI (fftw_mpi_) function, and before importing any wisdom for MPI problems, you must call:

     void fftw_mpi_init(void);

If FFTW threads support is used, however, fftw_mpi_init should be called after fftw_init_threads (see Combining MPI and Threads).  Calling fftw_mpi_init additional times (before fftw_mpi_cleanup) has no effect.

If you want to deallocate all persistent data and reset FFTW to the pristine state it was in when you started your program, you can call:

     void fftw_mpi_cleanup(void);

(This calls fftw_cleanup, so you need not call the serial cleanup routine too, although it is safe to do so.)  After calling fftw_mpi_cleanup, all existing plans become undefined, and you should not attempt to execute or destroy them.  You must call fftw_mpi_init again after fftw_mpi_cleanup if you want to resume using the MPI FFTW routines.

Using MPI Plans
---------------

Once an MPI plan is created, you can execute and destroy it using fftw_execute, fftw_destroy_plan, and the other functions in the serial interface that operate on generic plans (see Using Plans).

The fftw_execute and fftw_destroy_plan functions, applied to MPI plans, are collective calls: they must be called for all processes in the communicator that was used to create the plan.

You must not use the serial new-array plan-execution functions fftw_execute_dft and so on (see New-array Execute Functions) with MPI plans.  Such functions are specialized to the problem type, and there are specific new-array execute functions for MPI plans:

     void fftw_mpi_execute_dft(fftw_plan p, fftw_complex *in, fftw_complex *out);
     void fftw_mpi_execute_dft_r2c(fftw_plan p, double *in, fftw_complex *out);
     void fftw_mpi_execute_dft_c2r(fftw_plan p, fftw_complex *in, double *out);
     void fftw_mpi_execute_r2r(fftw_plan p, double *in, double *out);

These functions have the same restrictions as those of the serial new-array execute functions.  They are always safe to apply to the same in and out arrays that were used to create the plan.  They can only be applied to new arrarys if those arrays have the same types, dimensions, in-placeness, and alignment as the original arrays, where the best way to ensure the same alignment is to use FFTW's fftw_malloc and related allocation functions for all arrays (see Memory Allocation).  Note that distributed transposes (see FFTW MPI Transposes) use fftw_mpi_execute_r2r, since they count as rank-zero r2r plans from FFTW's perspective.

MPI Data Distribution Functions
-------------------------------

As described above (see MPI Data Distribution), in order to allocate your arrays, before creating a plan, you must first call one of the following routines to determine the required allocation size and the portion of the array locally stored on a given process.  The MPI_Comm communicator passed here must be equivalent to the communicator used below for plan creation.

The basic interface for multidimensional transforms consists of the functions:

     ptrdiff_t fftw_mpi_local_size_2d(ptrdiff_t n0, ptrdiff_t n1, MPI_Comm comm,
                                      ptrdiff_t *local_n0, ptrdiff_t *local_0_start);
     ptrdiff_t fftw_mpi_local_size_3d(ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t n2,
                                      MPI_Comm comm,
                                      ptrdiff_t *local_n0, ptrdiff_t *local_0_start);
     ptrdiff_t fftw_mpi_local_size(int rnk, const ptrdiff_t *n, MPI_Comm comm,
                                   ptrdiff_t *local_n0, ptrdiff_t *local_0_start);

     ptrdiff_t fftw_mpi_local_size_2d_transposed(ptrdiff_t n0, ptrdiff_t n1, MPI_Comm comm,
                                                 ptrdiff_t *local_n0, ptrdiff_t *local_0_start,
                                                 ptrdiff_t *local_n1, ptrdiff_t *local_1_start);
     ptrdiff_t fftw_mpi_local_size_3d_transposed(ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t n2,
                                                 MPI_Comm comm,
                                                 ptrdiff_t *local_n0, ptrdiff_t *local_0_start,
                                                 ptrdiff_t *local_n1, ptrdiff_t *local_1_start);
     ptrdiff_t fftw_mpi_local_size_transposed(int rnk, const ptrdiff_t *n, MPI_Comm comm,
                                              ptrdiff_t *local_n0, ptrdiff_t *local_0_start,
                                              ptrdiff_t *local_n1, ptrdiff_t *local_1_start);

These functions return the number of elements to allocate (complex numbers for DFT/r2c/c2r plans, real numbers for r2r plans), whereas the local_n0 and local_0_start return the portion (local_0_start to local_0_start + local_n0 - 1) of the first dimension of an ndims array that is stored on the local process.  See Basic and advanced distribution interfaces.  For FFTW_MPI_TRANSPOSED_OUT plans, the _transposed variants are useful in order to also return the local portion of the first dimension in the ndimstrans transposed output.  See Transposed distributions.  The advanced interface for multidimensional transforms is:

     ptrdiff_t fftw_mpi_local_size_many(int rnk, const ptrdiff_t *n, ptrdiff_t howmany,
                                        ptrdiff_t block0, MPI_Comm comm,
                                        ptrdiff_t *local_n0, ptrdiff_t *local_0_start);
     ptrdiff_t fftw_mpi_local_size_many_transposed(int rnk, const ptrdiff_t *n, ptrdiff_t howmany,
                                                   ptrdiff_t block0, ptrdiff_t block1, MPI_Comm comm,
                                                   ptrdiff_t *local_n0, ptrdiff_t *local_0_start,
                                                   ptrdiff_t *local_n1, ptrdiff_t *local_1_start);

These differ from the basic interface in only two ways.  First, they allow you to specify block sizes block0 and block1 (the latter for the transposed output); you can pass FFTW_MPI_DEFAULT_BLOCK to use FFTW's default block size as in the basic interface.  Second, you can pass a howmany parameter, corresponding to the advanced planning interface below: this is for transforms of contiguous howmany-tuples of numbers (howmany = 1 in the basic interface).

The corresponding basic and advanced routines for one-dimensional transforms (currently only complex DFTs) are:

     ptrdiff_t fftw_mpi_local_size_1d(
                  ptrdiff_t n0, MPI_Comm comm, int sign, unsigned flags,
                  ptrdiff_t *local_ni, ptrdiff_t *local_i_start,
                  ptrdiff_t *local_no, ptrdiff_t *local_o_start);
     ptrdiff_t fftw_mpi_local_size_many_1d(
                  ptrdiff_t n0, ptrdiff_t howmany,
                  MPI_Comm comm, int sign, unsigned flags,
                  ptrdiff_t *local_ni, ptrdiff_t *local_i_start,
                  ptrdiff_t *local_no, ptrdiff_t *local_o_start);

As above, the return value is the number of elements to allocate (complex numbers, for complex DFTs).  The local_ni and local_i_start arguments return the portion (local_i_start to local_i_start + local_ni - 1) of the 1d array that is stored on this process for the transform input, and local_no and local_o_start are the corresponding quantities for the input.  The sign (FFTW_FORWARD or FFTW_BACKWARD) and flags must match the arguments passed when creating a plan.  Although the inputs and outputs have different data distributions in general, it is guaranteed that the output data distribution of an FFTW_FORWARD plan will match the input data distribution of an FFTW_BACKWARD plan and vice versa; similarly for the FFTW_MPI_SCRAMBLED_OUT and FFTW_MPI_SCRAMBLED_IN flags. See One-dimensional distributions.

MPI Plan Creation
-----------------

Complex-data MPI DFTs
.....................

Plans for complex-data DFTs (see 2d MPI example) are created by:

     fftw_plan fftw_mpi_plan_dft_1d(ptrdiff_t n0, fftw_complex *in, fftw_complex *out,
                                    MPI_Comm comm, int sign, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_2d(ptrdiff_t n0, ptrdiff_t n1,
                                    fftw_complex *in, fftw_complex *out,
                                    MPI_Comm comm, int sign, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_3d(ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t n2,
                                    fftw_complex *in, fftw_complex *out,
                                    MPI_Comm comm, int sign, unsigned flags);
     fftw_plan fftw_mpi_plan_dft(int rnk, const ptrdiff_t *n,
                                 fftw_complex *in, fftw_complex *out,
                                 MPI_Comm comm, int sign, unsigned flags);
     fftw_plan fftw_mpi_plan_many_dft(int rnk, const ptrdiff_t *n,
                                      ptrdiff_t howmany, ptrdiff_t block, ptrdiff_t tblock,
                                      fftw_complex *in, fftw_complex *out,
                                      MPI_Comm comm, int sign, unsigned flags);

These are similar to their serial counterparts (see Complex DFTs) in specifying the dimensions, sign, and flags of the transform.  The comm argument gives an MPI communicator that specifies the set of processes to participate in the transform; plan creation is a collective function that must be called for all processes in the communicator.  The in and out pointers refer only to a portion of the overall transform data (see MPI Data Distribution) as specified by the local_size functions in the previous section.  Unless flags contains FFTW_ESTIMATE, these arrays are overwritten during plan creation as for the serial interface.  For multi-dimensional transforms, any dimensions > 1 are supported; for one-dimensional transforms, only composite (non-prime) n0 are currently supported (unlike the serial FFTW).  Requesting an unsupported transform size will yield a NULL plan.  (As in the serial interface, highly composite sizes generally yield the best performance.)

The advanced-interface fftw_mpi_plan_many_dft additionally allows you to specify the block sizes for the first dimension (block) of the ndims input data and the first dimension (tblock) of the ndimstrans transposed data (at intermediate steps of the transform, and for the output if FFTW_TRANSPOSED_OUT is specified in flags).  These must be the same block sizes as were passed to the corresponding local_size function; you can pass FFTW_MPI_DEFAULT_BLOCK to use FFTW's default block size as in the basic interface.  Also, the howmany parameter specifies that the transform is of contiguous howmany-tuples rather than individual complex numbers; this corresponds to the same parameter in the serial advanced interface (see Advanced Complex DFTs) with stride = howmany and dist = 1.

MPI flags
.........

The flags can be any of those for the serial FFTW (see Planner Flags), and in addition may include one or more of the following MPI-specific flags, which improve performance at the cost of changing the output or input data formats.

  * FFTW_MPI_SCRAMBLED_OUT, FFTW_MPI_SCRAMBLED_IN: valid for 1d transforms only, these flags indicate that the output/input of the transform are in an undocumented ``scrambled'' order.  A forward FFTW_MPI_SCRAMBLED_OUT transform can be inverted by a backward FFTW_MPI_SCRAMBLED_IN (times the usual 1/N normalization). See One-dimensional distributions.

  * FFTW_MPI_TRANSPOSED_OUT, FFTW_MPI_TRANSPOSED_IN: valid for multidimensional (rnk > 1) transforms only, these flags specify that the output or input of an ndims transform is transposed to ndimstrans.  See Transposed distributions.

Real-data MPI DFTs
..................

Plans for real-input/output (r2c/c2r) DFTs (see Multi-dimensional MPI DFTs of Real Data) are created by:

     fftw_plan fftw_mpi_plan_dft_r2c_2d(ptrdiff_t n0, ptrdiff_t n1,
                                        double *in, fftw_complex *out,
                                        MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_r2c_2d(ptrdiff_t n0, ptrdiff_t n1,
                                        double *in, fftw_complex *out,
                                        MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_r2c_3d(ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t n2,
                                        double *in, fftw_complex *out,
                                        MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_r2c(int rnk, const ptrdiff_t *n,
                                     double *in, fftw_complex *out,
                                     MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_c2r_2d(ptrdiff_t n0, ptrdiff_t n1,
                                        fftw_complex *in, double *out,
                                        MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_c2r_2d(ptrdiff_t n0, ptrdiff_t n1,
                                        fftw_complex *in, double *out,
                                        MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_c2r_3d(ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t n2,
                                        fftw_complex *in, double *out,
                                        MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_dft_c2r(int rnk, const ptrdiff_t *n,
                                     fftw_complex *in, double *out,
                                     MPI_Comm comm, unsigned flags);

Similar to the serial interface (see Real-data DFTs), these transform logically ndims real data to/from ndimshalf complex data, representing the non-redundant half of the conjugate-symmetry output of a real-input DFT (see Multi-dimensional Transforms). However, the real array must be stored within a padded ndimspad array (much like the in-place serial r2c transforms, but here for out-of-place transforms as well). Currently, only multi-dimensional (rnk > 1) r2c/c2r transforms are supported (requesting a plan for rnk = 1 will yield NULL).  As explained above (see Multi-dimensional MPI DFTs of Real Data), the data distribution of both the real and complex arrays is given by the local_size function called for the dimensions of the complex array.  Similar to the other planning functions, the input and output arrays are overwritten when the plan is created except in FFTW_ESTIMATE mode.

As for the complex DFTs above, there is an advance interface that allows you to manually specify block sizes and to transform contiguous howmany-tuples of real/complex numbers:

     fftw_plan fftw_mpi_plan_many_dft_r2c
                   (int rnk, const ptrdiff_t *n, ptrdiff_t howmany,
                    ptrdiff_t iblock, ptrdiff_t oblock,
                    double *in, fftw_complex *out,
                    MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_many_dft_c2r
                   (int rnk, const ptrdiff_t *n, ptrdiff_t howmany,
                    ptrdiff_t iblock, ptrdiff_t oblock,
                    fftw_complex *in, double *out,
                    MPI_Comm comm, unsigned flags);

MPI r2r transforms
..................

There are corresponding plan-creation routines for r2r transforms (see More DFTs of Real Data), currently supporting multidimensional (rnk > 1) transforms only (rnk = 1 will yield a NULL plan):

     fftw_plan fftw_mpi_plan_r2r_2d(ptrdiff_t n0, ptrdiff_t n1,
                                    double *in, double *out,
                                    MPI_Comm comm,
                                    fftw_r2r_kind kind0, fftw_r2r_kind kind1,
                                    unsigned flags);
     fftw_plan fftw_mpi_plan_r2r_3d(ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t n2,
                                    double *in, double *out,
                                    MPI_Comm comm,
                                    fftw_r2r_kind kind0, fftw_r2r_kind kind1, fftw_r2r_kind kind2,
                                    unsigned flags);
     fftw_plan fftw_mpi_plan_r2r(int rnk, const ptrdiff_t *n,
                                 double *in, double *out,
                                 MPI_Comm comm, const fftw_r2r_kind *kind,
                                 unsigned flags);
     fftw_plan fftw_mpi_plan_many_r2r(int rnk, const ptrdiff_t *n,
                                      ptrdiff_t iblock, ptrdiff_t oblock,
                                      double *in, double *out,
                                      MPI_Comm comm, const fftw_r2r_kind *kind,
                                      unsigned flags);

The parameters are much the same as for the complex DFTs above, except that the arrays are of real numbers (and hence the outputs of the local_size data-distribution functions should be interpreted as counts of real rather than complex numbers).  Also, the kind parameters specify the r2r kinds along each dimension as for the serial interface (see Real-to-Real Transform Kinds).  See Other Multi-dimensional Real-data MPI Transforms.

MPI transposition
.................

FFTW also provides routines to plan a transpose of a distributed n0 by n1 array of real numbers, or an array of howmany-tuples of real numbers with specified block sizes (see FFTW MPI Transposes):

     fftw_plan fftw_mpi_plan_transpose(ptrdiff_t n0, ptrdiff_t n1,
                                       double *in, double *out,
                                       MPI_Comm comm, unsigned flags);
     fftw_plan fftw_mpi_plan_many_transpose
                     (ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t howmany,
                      ptrdiff_t block0, ptrdiff_t block1,
                      double *in, double *out, MPI_Comm comm, unsigned flags);

These plans are used with the fftw_mpi_execute_r2r new-array execute function (see Using MPI Plans), since they count as (rank zero) r2r plans from FFTW's perspective.

MPI Wisdom Communication
------------------------

To facilitate synchronizing wisdom among the different MPI processes, we provide two functions:

     void fftw_mpi_gather_wisdom(MPI_Comm comm);
     void fftw_mpi_broadcast_wisdom(MPI_Comm comm);

The fftw_mpi_gather_wisdom function gathers all wisdom in the given communicator comm to the process of rank 0 in the communicator: that process obtains the union of all wisdom on all the processes.  As a side effect, some other processes will gain additional wisdom from other processes, but only process 0 will gain the complete union.

The fftw_mpi_broadcast_wisdom does the reverse: it exports wisdom from process 0 in comm to all other processes in the communicator, replacing any wisdom they currently have.

See FFTW MPI Wisdom.

FFTW MPI Fortran Interface
==========================

The FFTW MPI interface is callable from modern Fortran compilers supporting the Fortran 2003 iso_c_binding standard for calling C functions.  As described in Calling FFTW from Modern Fortran, this means that you can directly call FFTW's C interface from Fortran with only minor changes in syntax.  There are, however, a few things specific to the MPI interface to keep in mind:

  * Instead of including fftw3.f03 as in Overview of Fortran interface, you should include 'fftw3-mpi.f03' (after use, intrinsic :: iso_c_binding as before).  The fftw3-mpi.f03 file includes fftw3.f03, so you should not include them both yourself.  (You will also want to include the MPI header file, usually via include 'mpif.h' or similar, although though this is not needed by fftw3-mpi.f03 per se.)  (To use the fftwl_ long double extended-precision routines in supporting compilers, you should include fftw3f-mpi.f03 in addition to fftw3-mpi.f03. See Extended and quadruple precision in Fortran.)

  * Because of the different storage conventions between C and Fortran, you reverse the order of your array dimensions when passing them to FFTW (see Reversing array dimensions).  This is merely a difference in notation and incurs no performance overhead.  However, it means that, whereas in C the first dimension is distributed, in Fortran the last dimension of your array is distributed.

  * In Fortran, communicators are stored as integer types; there is no MPI_Comm type, nor is there any way to access a C MPI_Comm.  Fortunately, this is taken care of for you by the FFTW Fortran interface: whenever the C interface expects an MPI_Comm type, you should pass the Fortran communicator as an integer. (Technically, this is because you aren't actually calling the C functions directly. You are calling wrapper functions that translate the communicator with MPI_Comm_f2c before calling the ordinary C interface.  This is all done transparently, however, since the fftw3-mpi.f03 interface file renames the wrappers so that they are called in Fortran with the same names as the C interface functions.)

  * Because you need to call the local_size function to find out how much space to allocate, and this may be larger than the local portion of the array (see MPI Data Distribution), you should always allocate your arrays dynamically using FFTW's allocation routines as described in Allocating aligned memory in Fortran. (Coincidentally, this also provides the best performance by guaranteeding proper data alignment.)

  * Because all sizes in the MPI FFTW interface are declared as ptrdiff_t in C, you should use integer(C_INTPTR_T) in Fortran (see FFTW Fortran type reference).

  * In Fortran, because of the language semantics, we generally recommend using the new-array execute functions for all plans, even in the common case where you are executing the plan on the same arrays for which the plan was created (see Plan execution in Fortran). However, note that in the MPI interface these functions are changed: fftw_execute_dft becomes fftw_mpi_execute_dft, etcetera. See Using MPI Plans.

For example, here is a Fortran code snippet to perform a distributed L,M complex DFT in-place.  (This assumes you have already initialized MPI with MPI_init and have also performed call fftw_mpi_init.)

       use, intrinsic :: iso_c_binding
       include 'fftw3-mpi.f03'
       integer(C_INTPTR_T), parameter :: L = ...
       integer(C_INTPTR_T), parameter :: M = ...
       type(C_PTR) :: plan, cdata
       complex(C_DOUBLE_COMPLEX), pointer :: data(:,:)
       integer(C_INTPTR_T) :: i, j, alloc_local, local_M, local_j_offset

     !   get local data size and allocate (note dimension reversal)
       alloc_local = fftw_mpi_local_size_2d(M, L, MPI_COMM_WORLD, &
                                            local_M, local_j_offset)
       cdata = fftw_alloc_complex(alloc_local)
       call c_f_pointer(cdata, data, [L,local_M])

     !   create MPI plan for in-place forward DFT (note dimension reversal)
       plan = fftw_mpi_plan_dft_2d(M, L, data, data, MPI_COMM_WORLD, &
                                   FFTW_FORWARD, FFTW_MEASURE)

     ! initialize data to some function my_function(i,j)
       do j = 1, local_M
         do i = 1, L
           data(i, j) = my_function(i, j + local_j_offset)
         end do
       end do

     ! compute transform (as many times as desired)
       call fftw_mpi_execute_dft(plan, data, data)

       call fftw_destroy_plan(plan)
       call fftw_free(cdata)

Note that when we called fftw_mpi_local_size_2d and fftw_mpi_plan_dft_2d with the dimensions in reversed order, since a L,M Fortran array is viewed by FFTW in C as a M, L array.  This means that the array was distributed over the M dimension, the local portion of which is a L,local_M array in Fortran.  (You must not use an allocate statement to allocate an L,local_M array, however; you must allocate alloc_local complex numbers, which may be greater than L * local_M, in order to reserve space for intermediate steps of the transform.)  Finally, we mention that because C's array indices are zero-based, the local_j_offset argument can conveniently be interpreted as an offset in the 1-based j index (rather than as a starting index as in C).

If instead you had used the ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT) flag, the output of the transform would be a transposed M,local_L array, associated with the same cdata allocation (since the transform is in-place), and which you could declare with:

       complex(C_DOUBLE_COMPLEX), pointer :: tdata(:,:)
       ...
       call c_f_pointer(cdata, tdata, [M,local_L])

where local_L would have been obtained by changing the fftw_mpi_local_size_2d call to:

       alloc_local = fftw_mpi_local_size_2d_transposed(M, L, MPI_COMM_WORLD, &
                                local_M, local_j_offset, local_L, local_i_offset)