Complex arrays

The story for arrays of complex numbers is, well, a bit more complex. There are two possible aways of arranging complex numbers (the indexing below is for data centering):

Interleaved: In this approach, the real and imaginary parts are stored in successive array locations in arrays of [$x$,$y$,$z$]:

Re[0,0,0], Im[0,0,0], Re[1,0,0], Im[1,0,0], $\ldots$

This is the most common approach.

Split: In this approach, the entire real array part is stored contiguously, followed by the entire imaginary part in arrays of [$x$,$y$,$z$]:

Re[0,0,0], Re[1,0,0], Re[2,0,0], $\ldots$,
Im[0,0,0], Im[1,0,0], Im[2,0,0], $\ldots$

This approach turns out to be slightly faster on some parallel computers (such as the dist_fft routines on an Apple cluster; see Sec. 6.2).

To make code that will work in either of these two conventions, you must follow some special mechanisms to work with complex arrays:

• The type dm_array_index_t is presently set to a 32 bit unsigned integer which works for arrays up to $1024^{3}$, but by using dm_array_index_t you allow for code to be modified to work with larger arrays and 64 bit indexing in the future by a simple re-definition of the typedef. Also, note that sizeof(dm_array_complex) will not necessarily give you the proper array size for storage because for split data it gives you the length of two memory references; see Sec. 1.4.3 on how to allocate memory for complex arrays.

• To manipulate a pixel in the complex array, do not refer to it as
  real_part = *(complex_array+2*(ix+iy*nx+iz*nx*ny));
imaginary_part = *(complex_array+2*(ix+iy*nx+iz*nx*ny)+1);
  but instead use the macros c_re and c_im for access:
  real_part = c_re(complex_array,(ix+iy*nx+iz*nx*ny));
imaginary_part = c_im(complex_array,(ix+iy*nx+iz*nx*ny));

Also note that the typedefs of dm_array_real and dm_array_complex change between float and double according to the presence or absence of the #define variable DM_ARRAY_DOUBLE.