Question on F-order PlanNd

[leofang]

@grlee77 I have a quick question on lines like this:

cupy/cupy/fft/fft.py

Lines 287 to 288 in fa83171

	if order == 'F':
	plan_dimensions = plan_dimensions[::-1]

I am having trouble understanding why it works. Could you explain to me? Thanks!

[grlee77]

If you look at the notes section for PlanNd, you will see the following (copied from the CUDA docs):

    in 3D:
    input[b * idist + ((x * inembed[1] + y) * inembed[2] + z) * istride]
    output[b * odist + ((x * onembed[1] + y) * onembed[2] + z) * ostride]

Here:
x = axis 0
y = axis 1
z = axis 2
Axis "z" has the smallest stride, so the indexing described there is assuming order='C' (elements along the last axis are adjacent in memory). There is no option in CUFFT to specify "order=F", but if will give the reversed shape, the indexing will come out correct. This is a little tricky to think through, but that is the basic idea.

A concrete example to illustrate this:

shape = (32, 16, 8)
print(np.ones(shape, order='C').strides)
(1024, 64, 8)

print(np.ones(shape, order='F').strides)
(8, 256, 4096)

print(np.ones(shape[::-1], order='F').strides)
(8, 64, 1024)

Note that order 'F' with reversed shape at least has the same numerical values for the strides (although the order is reversed).

[grlee77]

I think someone had recommended this method of reversing the shape somewhere on the NVIDIA forums years ago, but I don't recall exactly where. A quick google search found the following section of a book recommending swapping nx, ny for 2D FFTs, although that is not where I originally saw the suggestion.

[grlee77]

Actually, it looks like the shape reversal used to be mentioned in the CUFFT docs, but I don't see it in the current ones. In this older version it states:

"For higher‐dimensional transforms (2D and 3D), CUFFT performs FFTs in row‐major or C order. For example, if the user requests a 3D transform plan for sizes X, Y, and Z, CUFFT transforms along Z, Y,
and then X. The user can configure column‐major FFTs by simply changing the order of size parameters to the plan creation API functions.
"

[leofang]

Forgot to reply here: Thanks a lot, Gregory! I was able to see that when the indices [z][y][x] is reversed to [x][y][z] for F-order arrays, the flatted index [b * idist + ((x * inembed[1] + y) * inembed[2] + z) * istride] would be 0, 1, 2, 3, ... as the original 3D indexes increase indeed. But I wasn't so sure if this alone is enough for correctness.

Now, I am suspecting I need to think the indexing harder for R2C/C2R. In my PR I just found F-order inputs sometimes are transformed wrong. It may not be as simple as in C2C...?

[leofang]

Leaving a note for myself: Still trying to figure out how reversing array dimensions works for R2C & C2R transforms of F-order arrays. I thought it should be just a straightforward reverse, but doesn't seem like it.

I had a sense that cuFFT works quite similarly as FFTW, so I looked up the FFTW doc and found the same suggestion for dimension reverse:
http://www.fftw.org/fftw3_doc/Reversing-array-dimensions.html#Reversing-array-dimensions
So I have no clue now and will need to revisit later...

[leofang]

I think cuFFT R2C/C2R is probably buggy for 2D transform + F-order...For 3D the shape reversal trick works just fine. 😞

[leofang]

Found something interesting. Let me redirect the discussion back to my PR cupy#3102.

[leofang]

For my own note: I found this on FFTW's doc:

We should remind the user that the separable product of 1d transforms along each dimension, as computed by FFTW, is not always the same thing as the usual multi-dimensional transform.

See http://www.fftw.org/fftw3_doc/Multi_002ddimensional-Transforms.html#Multi_002ddimensional-Transforms.

"Separable product of 1d transforms along each dimension" is what's done in NumPy. So I wonder if this is the real cause of mismatch between NumPy and cuFFT (suppose cuFFT follows closely what's done in FFTW)...