Re: [eigen] A complex FFT for Eigen

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

Arf good job! One typo, it's ei_pstore instead of ei_store, and
ei_pmadd instead of ei_madd (just in case Tim tries this code).

One thing though, is it premature optimization to write it in such a
lowlevel way even before we have global vectorization support for
complex numbers... I mean, your code is perfect for complex<float> but
leaves complex<double> outside in the cold. Maybe it'd be easier to
get a solution that works for both once we have global vectorization
support for complex. Maybe the first thing to do would be to add
specialisations of PacketMath.h for complex types, so we don't even
need to have our own complex class.

Cheers,
Benoit

2008/12/2 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>
> all of this looks great,
>
> On Mon, Dec 1, 2008 at 11:31 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
> wrote:
>>
>> 2008/12/1 Tim Molteno <tim@xxxxxxxxxxxxxxxxxxx>:
>>
>>                for (unsigned i=0; i<N2; i++)
>>                {
>>                        std::complex<T> temp(data[i+N2]*w);
>>                        data[i+N2] = data[i]-temp;
>>                        data[i] += temp;
>>                        w += w*wp;
>>                }
>>
>> And I don't see any way to use Eigen expressions here, or to vectorize
>> that loop as such (again, no contiguous access).
>
> what about the following (let's assume the size of a packet is 2):
>
> const int PacketSize = 2;
> Scalar _w[PacketSize];
> _w[0] = 0;
> for (int k=1; k<PacketSize; ++k)
>   _w[k] = _w[k-1] + _w[k-1] * wp;
> Packet w = ei_pload(_w);
> Packet wp2 = ei_pset1(wp);
> for (int i=0; i<N2; i+=PacketSize)
> {
>   Packet tmp = ei_pmul(ei_pload(data+i+N2), w);
>   Packet di = ei_pload(data+i);
>   ei_store(data+i+N2, ei_psub(di, tmp));
>   ei_store(data+i,       ei_padd(di, tmp));
>   w = ei_madd(w, wp2, w);
> }
>
> I think this should work whatever PacketSize is (including the case when the
> vectorization is disabled, i.e., PacketSize = 1)
>
> gael.
>

---