FFT produces unexpected transposed result when either input or output is Fortran-ordered
cmey opened this issue · 1 comments
I built a small test code to reproduce the issue: https://gist.github.com/cmey/e45ed9dc2b04bc216dcf032d50aeae99
It fails when either the input or the output array is Fortran-ordered, example output:
reikna: array([[ 1.4e+00+0.j, 8.3e-01+0.j],
[ -6.1e-01+0.j, -1.1e-03+0.j]], dtype=complex64)
np_fft: array([[ 1.4e+00+0.j, -6.1e-01+0.j],
[ 8.3e-01+0.j, -1.1e-03+0.j]])
Yes, currently reikna only supports C-ordered arrays. Of course, it should, at the very least, raise an exception if a Fortran-ordered array is passed.
Technically, all global memory loads/stores are done through load_combined_idx() and store_combined_idx(), so the Fortran-order support is easy to add there. The problem is that many computations, including FFT, rely on the last dimension being the innermost, so they can take advantage of memory coalescing. Another option is to recognize F-ordered arrays at the point of creation of Type objects and reverse the shape if necessary, so that computations work with C-ordered shapes. I am not sure what is the best way to approach this, I'll have to think about it.