[Gsll-devel] Possible use of user-supplied vector sizes and offsets for FFT

[Sumant Oemrawsingh]

Hi Liam,

As I briefly indicated on #lisp and contrary to our earlier belief, there _is_
at least one case in which a user would like to supply the stride, vector size
n and offset to the FFT functions. Namely, this is for the 2-dimensional (more
generally of course, N-dimensional) Fourier transform.

GSL only provides one-dimensional FFTs. Naturally, an N-dimensional FFT can be
easily implemented from a one-dimensional function. From saying this, you can
probably already immediately see what I mean, but nevertheless, to be sure
there's no mistake here, I'll give a more concrete example.

Let us approach the problem in C, and consider a 2D array of size a x b, where
a is the number of rows and b is the number of columns. This is a linear,
contiguous array in C (a "raw" C array). A 2D FFT can be performed by first
FFT'ing one dimension and then the other.

Thus, we first FFT each row, and then each column. The order actually doesn't
matter. In C, this would be handled by first performing a times the (forward)
FFT, thus applying an FFT to each row, where the vector size n = b, the stride
= 1 and each consecutive FFT call supplies the vector with an offset with a
multiple of b.

Next up is the FFT of each column. In C we perform b times the (forward) FFT,
thus applying an FFT to each column where the vector size n = a, the stride =
b, and each consecutive FFT call supplies the vector with an offset equal to
the column number.

For N-dimensional arrays, simply lather, rinse and repeat.

This shows the necessity of the optional (keyword?) argument that indicates
the size of the vector to be FFT'd, the stride (which is already in there) and
also the offset (which is missing at the moment). This also shows that the
Fourier transform environment that is in fft-interface.lisp is very nice to
have, when dealing with non-radix-2 array lengths. Each consecutive call when
FFT-ing in one direction only needs to allocate the workspace and wavetable
once and, more importantly, calculate the wavetable once, while the FFT is
repeated on each row (or column).

I have to admit that I haven't looked at how multi-dimensional marrays or more
correctly, marrays representing tensors of rank > 1, from the GSLL point of
view. As I understand, the interface allows for a truly multi-rank approach
(using maref with multiple indices), as with the normal CL arrays and aref.

In my opinion, a nice approach would be to make the whole FFT by default work
on the entire marray, thus performing multi-dimensional FFTs when marrays have
higher ranks. That way, we can hide the messing about with strides, sizes and
offsets from the user. Just think, loading an image into a marray, and just
calling forward-fourier-transform, and instantly (depending on the speed of
your computer...) have a Fourier-space image!

On the other hand, while such an interface would be nice for the casual user
(and probably for most normal users) I don't know if there are "power users"
that (would want to use GSL's FFT functions, and) would require direct control
over such parameters.

As I may have mentioned before on #lisp, I have used GSL's FFT some time ago,
as well as FFTW briefly. FFTW has functions for higher-dimensional FFTs, and
deducing from that interface, I would say that an FFT function that "does the
right thing (TM)" depending on the rank of the marray would not be so bad.

Anyway, sorry for the long post. I hope my explanations/considerations/ideas
are useful and I would not mind implementing some/most of this in an
experimental interface, possibly with some help with respect to some of the
fancy internals of GSLL.

Please let me know what you think (and any suggestions from users of GSLL and
FTT, I'm looking at you, Mirko...) are greatly appreciated!

Thanks,
Sumant