Re: [eigen] A complex FFT for Eigen

["Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>]

Hi,

I agree with what you say.

I agree that it would be useful to have a FFT module that has both a
builtin implementation (Tim's code is the beginning of it) and that
allows to use FFTW as optional backend.

As Keir notes, FFTW however is GPL so at least GPL-incompatible users
would always be using our builtin implementation.

We need extension to dynamic-size; of course that can't be done with
template metaprogramming but would involve a dynamic loop; we could
still try to "peel" it partly into a fixed-size version. Gael's idea
is exactly that and of course it applies not only to dynamic size but
also to too large fixed sizes.

Regarding complex numbers vectorization: this means that we'd need to
provide our own ei_complex type... and that we'd need to become much
more careful than we are now as our current code, at many places, is
based around the assumption that PacketSize==1 implies no
vectorization; also some of our vectorized paths (like in Dot.h)
assume real numbers as for now we don't vectorize complex numbers.

Cheers,
Benoit

2008/11/28 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>
> Hi,
>
> here are my thoughts.
>
> First of all, I'm 100% OK to provide a Fourier module in Eigen. My issues
> are more about how to do it. Looking at FFTW benchmarks, it is pretty clear
> to me that we *must* provide support for it. So, in the same vein of the
> other modules, I would provide a default implementation based on Tim's
> lightweight source code such that we can make the dependency to FFTW
> optional for most use cases of this module. Of course this requires to
> extend a bit Tim's code to support at least: real data, dynamic sizes, 1D
> and 2D data (I guess all are trivial to add). IMO non pot sizes is not
> mandatory.
>
> About current Tim's code, what about adding a specialization of
> Radix2_Decimation<> for sizes greater than, let's say 32, which would use
> normal recursive calls untill the size is 32. This would automatically add
> support for dynamic sizes, and significantly reduce the size of the
> generated code.
>
> About the vectorization of complex types, this is indeed perfectly doable.
> However note that for std::complex<double> the size of the packet would be 1
> ! So, somehow, this means that for double it is rather the std::complex's
> operators that should be updated to use the SSE instruction set....
>
> Gael.
>
>
> On Fri, Nov 28, 2008 at 3:12 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
> wrote:
>>
>> 2008/11/27 Tim Molteno <tim@xxxxxxxxxxxxxxxxxxx>:
>> > I'd love to vectorize this better. I'll look into this, but ran into the
>> > issue of addressing sub-vectors in-place. I could not work out how to do
>> > it
>> > in Eigen. How do I carry out operations like
>> >
>> > data[1:N] += data[N+1:2N]
>>
>> assuming that data is a vector of length 2N :
>>
>> data.start<N>() += data.end<N>();
>>
>> otherwise:
>>
>> data.start<N>() += data.segment<N>(N);
>>
>> (I think that segment only got that name after beta1 so you should use
>> trunk and "make doc" instead of using the online documentation).
>>
>> >
>> > In other words, can I create a vector object that refers into part of
>> > another vector object?
>> >
>> > I tried using MatrixBlock objects, but could not get it to work?
>>
>> Yes Eigen::Block is exactly what you need here and the above-mentioned
>> start/end/segment methods return objects of type Block.
>>
>> > I think that FFT's are inherently non-local operations.
>>
>> OK, so be it, then :)
>>
>> > There might be
>> > better ways of doing this, but as it stands this code is pretty fast for
>> > an
>> > FFT and extremely small and lightweight, so I'd rather vectorize the
>> > apply
>> > function first as this is where most of the computation is happening.
>>
>> Good to know.
>> As I said for now Eigen doesn't vectorize operations involve complex
>> numbers but I don't see anything making it impossible in theory, I'll
>> think about that when I can.
>>
>> So, from Matthias's and your answers, I gather that there is no
>> generalization to non-complex numbers or to non-power-of-two sizes.
>> What about dynamic-size: are there, in your experience, use cases
>> where the size N isn't known at compile-time and it's be useful to
>> have a FFT able to handle any size N as a runtime parameter? If yes,
>> that's also something that we want.
>>
>> More nitpicking:
>> Why are fft() and ifft() methods of a class CFFT ? I'd suggest either
>> making them static methods, or global functions. In Eigen, we tend to
>> do the latter as it makes for shorter calling syntax.
>>
>> Let's give us a bit more time to determine whether we want FFT inside
>> Eigen -- don't want to sound un-enthusiastic, it's rather that I'm
>> currently busy and I don't want to make the wrong decision.
>>
>> Anybody reading this, opinions?
>>
>> Cheers,
>> Benoit
>>
>> >
>> > Cheers
>> >
>> > Tim
>> >
>> > ---
>> >
>> >
>>
>> ---
>>
>
>

---