Re: [eigen] A complex FFT for Eigen |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] A complex FFT for Eigen*From*: "Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>*Date*: Tue, 2 Dec 2008 14:10:33 +0100*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:received:message-id:date:from:to :subject:in-reply-to:mime-version:content-type :content-transfer-encoding:content-disposition:references; bh=MW/+f7wSPhqUv+auGJvqwQ97fEcF4c9gaw7iFC6BWeA=; b=ks3hpOUML9SZ5xugTSpOxtxJmQDMUKqRntXHJ077PWdkeMPyiFzhuniuear1MYC2P3 Avurf4jr5/S3UWTXimUL4ChKPLitx5q22TQLDzY5j2Kz+VoSFJhWYnnq6Se0R7me40iO RHeSMHYT9Ppwqa3UC/NKaNuHEfRuX2OZZi2QM=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=message-id:date:from:to:subject:in-reply-to:mime-version :content-type:content-transfer-encoding:content-disposition :references; b=BPluBJ9FQq0bt7Rzx95t72X076PHD8G93YQikLHSIMR6Qh6PxiVyFPuTl1f+L/9YtV QscTSAwBY5UWXY/NigcIKN7VtIxh6xPfBhQZ4OoG8YGXrSScvvukXCrw6B3X7dDDjz/n pOORy59hj/cSH8r7fERJxFK1xc9h7DyMc32Wc=

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. > ---

**Follow-Ups**:**Re: [eigen] A complex FFT for Eigen***From:*Gael Guennebaud

**References**:**Re: [eigen] A complex FFT for Eigen***From:*Benoit Jacob

**Re: [eigen] A complex FFT for Eigen***From:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] A complex FFT for Eigen** - Next by Date:
**Re: [eigen] A complex FFT for Eigen** - Previous by thread:
**Re: [eigen] A complex FFT for Eigen** - Next by thread:
**Re: [eigen] A complex FFT for Eigen**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |