Why do the strides work in this way?

[ceruttids]

I'm opening this issue in a bid to get some help on a vexing code problem in my project, which incorporates PocketFFT. If this is not the right place to pose a question, I'm happy to take it offline with someone who would offer to help.

The principal problem is that I have Fortran-ordered multi-dimensional arrays, (column-major, the contiguous memory lists all column elements consecutively, in three dimensions {x, y, z} the column being a "vertical" column against x / y "slabs") , whereas PocketFFT expects C order (row-major). To help frame the issue, I have been able to carry out forward and backward FFT "round trips" with such data and access the individual elements in an order that displays them in the same output format as GNU Octave performing fftn(x) on the same multi-dimensional data. (Whether Octave uses Fortran- or C-ordered arrays, I am entering the data in a script so the rows, slabs, and columns are all implied by the layout of data in the script.) My question is more about why what I'm doing works--I don't want to be caught getting the right answers for the wrong reasons.

What I do for a given real-to-half complex problem with { m, n, p } elements:

• Make the input shape_in of the problem a std::vector<size_t> of { p, n, m }
• Under this logic of reversing the order of lengths fed to the API, the output shape of the problem would seem to be { p, n, m/2 + 1 }. This is not part of the PocketFFT API calls, however.
• The input strides stride_in work out to { m x n x R, m x R, R }, where R is sizeof(your_real_type).
• The output strides stride_out, which are part of the actual API call, then work out as { (m / 2 + 1) x n x C, (m / 2 + 1) x C, C }, where C is sizeof(std::complex<your_real_type>).

If I then run it like so:

r2c(shape_in, stride_in, stride_out, { 0, 1, 2 }, FORWARD, your_real_signal_data, your_complex_xform_data, 1.0);
c2r(shape_in, stride_out, stride_in, { 0, 1, 2 }, BACKWARD, your_complex_xform_data, your_real_signal_data, 1.0 / (m * n * p));

The resulting signal will make the correct round-trip. I can also query the intermediate complex_xform_data to produce what looks like the correct result (matching GNU Octave running the same problem--I got over the roundoff problem by rounding all input elements to the nearest 1.0e-8 and printing the input signal as %12.8lf so I could copy it into the Octave script). It took me a lot of effort to solve the logical Rubik's cube.

The insight, it seems, was that I need to access the elements of the transformed signal as if the result were (or, because it is?) now a problem with shape { m / 2 + 1, n, p }, but again in Fortran order. This pertains to the result of just the forward FFT, not the round trip. (In order to do a convolution, I will need to get in there, between forward and reverse FFTs, and do some element-wise operations. That will require me to know what elements of the frequency space I'm accessing.) In order to access element { x, y, z } of the transformed result, with x varying in the range [0, m) and y between [0, n) and z between [0, p), here's the pseudocode:

getFrequencyData(const size_t x_pos, const size_t y_pos, const size_t z_pos) {

  // Apply periodic considerations.  If the FFT is real-to-complex, the following filter will       
  // detect a request for information that might be outside the explicit data array.  A             
  // complex-to-complex transform will have frequencies for each piece of real data.                
  size_t x_act = x_pos - (shape_in[0] * (x_pos / shape_in[0]));
  size_t y_act = y_pos - (shape_in[1] * (y_pos / shape_in[1]));
  size_t z_act = z_pos - (shape_in[2] * (z_pos / shape_in[2]));
  size_t xyz_act;
  if (x_act >= shape_out[2]) {
    x_act = shape_in[0] - x_act;
    if (y_act > 0) {
      y_act = shape_in[1] - y_act;
    }
    if (z_act > 0) {
      z_act = shape_in[2] - z_act;
    }
  }

  // Determine the array index to access.
  size_t xyz_act = (((z_act * shape_out[1]) + y_act) * shape_out[2]) + x_act;
  result = z_frequency[xyz_act];

  // Take the complex conjugate, if necessary.
  if (x_act >= shape_out[2]) {
    result.y = -result.y;
  }
  return result;
}

I'm hoping that someone is here who can confirm my reasoning. This has been quite a puzzle, as I said, but I think I finally have it. I'm finding that PocketFFT and cuFFT both follow the same data layouts, including for real-to-half complex transforms, so I've been able to extend my own code's API to subsume both packages for FFTs in C++ and CUDA, respectively. With a few more improvements I should be able to leverage the class containing code such as the above without having to go back through this degree of pain every time I'm faced with a new FFT. It will also internalize the allocation of out-of-place data, if required, and even automate repeated calls to PocketFFT to mimic the behavior of batched FFT calls in cuFFT. Very glad to have your code holding up the C++ side of things, providing a critical CPU-only run mode and the ability to write other parts of the software that perform in serial C++ what we want all of our CUDA kernels to be doing up to 2000x faster.

[mreineck]

I'm not absolutely sure I understand what the problem is, but perhaps I can bring a little more clarity by saying that pocketFFT does, in fact, not require C-ordered arrays but can directly work on Fortran-ordered arrays as well (or any other ordering, sub-array views etc.)
If you have a 3D array with (nx, ny, nz) entries in Fortran ordering, you simply supply pocketFFT with a shape vector of (nx,ny,nz), and a stride vector of (S, nx*S, nx*ny*S), where S is the size o your data type in bytes.
Conversely, if the array was in C ordering instead, the shape would stay the same, but the strides would be (ny*nz*S, nz*S, S). pocketFFT was designed with the goal to minimize the pain caused by different memory ordering schemes, because I wanted it to be easily usable from different languages.
Long story short: you don't have to convert your (n,m,p) to a (p,n,m) array; just continue to view it as is natural to the environment you are calling from, and only adjust the strides.
Of course, in your specific case (r2c/c2r) there is the additional complication of the input and output arrays having different sizes along one axis, which doesn't make things easier. This is why I personally perer the Hartley transform in this case ... but since you need to stay compatible with cufft, this won't help you :-/

[ceruttids]

Thanks, Martin. I understand that PocketFFT does not require C ordering. What surprises me about your post is that I did not know I could supply the same shape vector but reverse / invert the strides. This makes sense but all of the advice I had seen online (and from colleagues) was to reverse the shape vector and alter the strides as I had done above.

My question is more focused on how I access data in the transformed result from a c2c or r2c operation. Hence, the code I posted. I'm not re-arranging the data in the actual array. I'd like to arrange the shape and stride vectors in the most transparent way to perform FFTs on Fortran-ordered multi-dimensional arrays, but however those shape and stride vectors are calculated that will now be encapsulated in the FFTStage class I've made to manage work with either PocketFFT or cuFFT (or hipFFT and other packaged from alternate HPC libraries if I get there). However the strides work out, I will need to perform work on transformed result, and while the above code will get inlined and streamlined in any major kernel / function that does those operations I'm posting to see if anyone has commentary on how I'm doing that indexing to get element { x, y, z } of the frequencies for a problem of { m, n, p } elements.

[mreineck]

Thanks, I think I understand the question better now, and I'm sorry this is not documented more clearly directly in the source code!

The insight, it seems, was that I need to access the elements of the transformed signal as if the result were (or, because it is?) now a problem with shape { m / 2 + 1, n, p }, but again in Fortran order.

This is indeed the case. This is mentioned as the last bullet point under "General constraints on arguments" in README.md, but I admit that this is not the most intuitive place to look for once you have started digging into the code itself.
The code generally has too little function-level documentation, since I intended it mainly as a kind of black-box engine behind a Python interface in scipy; all the interface documentation is present in the docstrings there.

[mreineck]

The choice of output array dimensions for the r2c transform is explained in the second paragraph of https://www.fftw.org/fftw3_doc/Multi_002ddimensional-Transforms.html. As far as I know, all FFT libraries implementing r2c use this format, potentially with slight generalizations. (For example, in FFTW the halved axis is always the last axis of the array, in pocketfft it is the last axis to be transformed, which may be any axis.)

[ceruttids]

Before I head off for the weekend, just wanted to say thanks for your contributions here. This is very helpful. The correspondence between PocketFFT's supported functionality and that of cuFFT is strong, save for one little detail about the availability of in-place real-to-complex transforms. I'll craft my wrapper to prevent developers from taking a wrong step in the C++ layer. For the GPU FFTs, I found in late 2023 that the global cache size is truly performance limiting (as one might expect), but that the newer GPUs (Hopper, Lovelace, and later) with much larger global caches smooth things out a great deal, even for problems that technically spill out of the cache. That's another thing I'm designing my API to automate, the batching of multiple FFTs insofar as the GPU cache can handle it.

[mreineck]

I'm aware that the lack of in-place c2r and r2c is painful for many users, and the decision not to support it was a difficult one. The problem is that there is no strictly language conforming way to make it work if different types are used to represent real and complex data, since C and C++ both state "if two pointers point to different data types (except char ...), memory accesses through those pointers can be rearranged arbitrarily, because they cannot (by definition in the language standard) point to overlapping memory locations".
In practice things work out properly in practically all situations, though. Maybe I should add some kind of "expert cheat mode" for this, but I'm not sure.

[mreineck]

I looked at the sources again, and in-place r2c/c2r should actually work, under the following circumstances:

• the first elements of input and output arrays must have exactly the same address
• the strides along axis or the last entry of axes must be exactly the size of one array entry for both input and output arrays
• the strides along all other axes must be identical (in bytes)

This is basically the same set of requirements that FFTW imposes.

I tried this on the pocketfft successor (a part of https://github.com/mreineck/ducc) and it worked well. I didn't have time for systematic tests with pocketfft itself though.