source: Sophya/trunk/SophyaExt/IFFTW/fftw2server.cc@ 3123

Last change on this file since 3123 was 3000, checked in by ansari, 19 years ago

Passage a FFTW3 + suppression const des arguments FFTForward/Backward cmv+Reza 03/07/2006
File size: 11.0 KB
Line 
1#include "FFTW/fftw.h"
2#include "FFTW/rfftw.h"
3
4
5namespace SOPHYA {
6
7class FFTWServerPlan {
8public:
9 FFTWServerPlan(int n, fftw_direction dir, bool fgreal=false);
10 FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal=false);
11 ~FFTWServerPlan();
12 void Recreate(int n);
13 void Recreate(int nd, int * nxyz);
14
15 static void FillSizes(const BaseArray & in, int * sz);
16
17 int _n; // Array dimension for 1-d arrays
18 int _nd; // Nb of dimensions for n-d arrays
19 int _nxyz[MAXND_FFTW]; // Array dimensions for n-d arrays
20 fftw_direction _dir;
21
22 fftw_plan p;
23 rfftw_plan rp;
24 fftwnd_plan pnd;
25 rfftwnd_plan rpnd;
26
27};
28
29} // Fin du namespace
30
31FFTWServerPlan::FFTWServerPlan(int n, fftw_direction dir, bool fgreal)
32{
33 if (n < 1)
34 throw ParmError("FFTWServerPlan: Array size <= 0 !");
35 p = NULL;
36 rp = NULL;
37 pnd = NULL;
38 rpnd = NULL;
39 _nd = -1;
40 for(int k=0; k<MAXND_FFTW; k++) _nxyz[k] = -10;
41 _n = n;
42 _dir = dir;
43 if (fgreal) {
44 rp = rfftw_create_plan(n, dir, FFTW_ESTIMATE);
45 if (rp == NULL)
46 throw AllocationError("FFTWServerPlan: failed to create plan (rp) !");
47 }
48 else {
49 p = fftw_create_plan(n, dir, FFTW_ESTIMATE);
50 if (p == NULL)
51 throw AllocationError("FFTWServerPlan: failed to create plan (p) !");
52 }
53
54}
55
56FFTWServerPlan::FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal)
57{
58 int k;
59 if (nd > MAXND_FFTW)
60 throw ParmError("FFTWServerPlan: Array rank (nd) > MAXND_FFTW !");
61 p = NULL;
62 rp = NULL;
63 pnd = NULL;
64 rpnd = NULL;
65 _n = -10;
66 _nd = nd;
67 for(k=0; k<nd; k++) {
68 if (nxyz[k] < 1)
69 throw ParmError("FFTWServerPlan: One of the Array size <= 0 !");
70 _nxyz[k] = nxyz[k];
71 }
72 for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
73 _dir = dir;
74 if (fgreal) {
75 rpnd = rfftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
76 if (rpnd == NULL)
77 throw AllocationError("FFTWServerPlan: failed to create plan (rpnd) !");
78 }
79 else {
80 pnd = fftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
81 if (pnd == NULL)
82 throw AllocationError("FFTWServerPlan: failed to create plan (pnd) !");
83 }
84
85}
86
87FFTWServerPlan::~FFTWServerPlan()
88{
89 if (p) fftw_destroy_plan(p);
90 if (rp) rfftw_destroy_plan(rp);
91 if (pnd) fftwnd_destroy_plan(pnd);
92 if (rpnd) rfftwnd_destroy_plan(rpnd);
93}
94
95void
96FFTWServerPlan::Recreate(int n)
97{
98 if (n < 1)
99 throw ParmError("FFTWServerPlan::Recreate(n) n < 0 !");
100 if (_nd > 0)
101 throw ParmError("FFTWServerPlan::Recreate(n) Multi-dimensional plan ! > 0 !");
102 if (n == _n) return;
103 _n = n;
104 if (p) {
105 fftw_destroy_plan(p);
106 p = fftw_create_plan(n, _dir, FFTW_ESTIMATE);
107 if (p == NULL)
108 throw AllocationError("FFTWServerPlan::Recreate failed to create plan (p) !");
109 }
110 else {
111 rfftw_destroy_plan(rp);
112 rp = rfftw_create_plan(n, _dir, FFTW_ESTIMATE);
113 if (rp == NULL)
114 throw AllocationError("FFTWServerPlan::Recreate failed to create plan (rp) !");
115 }
116
117}
118
119void
120FFTWServerPlan::Recreate(int nd, int * nxyz)
121{
122 if (_n > 0)
123 throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) 1-dimensional plan !");
124 int k;
125 if (nd == _nd) {
126 bool samepl = true;
127 for (int k=0; k<nd; k++)
128 if (nxyz[k] != _nxyz[k]) samepl = false;
129 if (samepl) return;
130 }
131 if (nd > MAXND_FFTW)
132 throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) Array rank (nd) > MAXND_FFTW !");
133 _nd = nd;
134 for(k=0; k<nd; k++) {
135 if (nxyz[k] < 1)
136 throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) One of the Array size <= 0 !");
137 _nxyz[k] = nxyz[k];
138 }
139 for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
140 if (pnd) {
141 fftwnd_destroy_plan(pnd);
142 pnd = fftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
143 if (pnd == NULL)
144 throw AllocationError("FFTWServerPlan::Recreate failed to create plan (pnd) !");
145 }
146 else {
147 rfftwnd_destroy_plan(rpnd);
148 rpnd = rfftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
149 if (rpnd == NULL)
150 throw AllocationError("FFTWServerPlan::Recreate failed to create plan (rpnd) !");
151 }
152
153}
154
155void
156FFTWServerPlan::FillSizes(const BaseArray & in, int * sz)
157{
158 int k1 = 0;
159 int k2 = 0;
160 for(k1=in.NbDimensions()-1; k1>=0; k1--) {
161 sz[k2] = in.Size(k1); k2++;
162 }
163}
164
165/* --Methode-- */
166//! Constructor - will always preserve input with FFTW2
167FFTWServer::FFTWServer(bool /*preserve_input*/)
168 : FFTServerInterface("FFTServer using FFTW package")
169 , ckR8("FFTWServer - ", true, false) ,
170 _preserve_input(true)
171
172{
173 _p1df = NULL;
174 _p1db = NULL;
175 _pndf = NULL;
176 _pndb = NULL;
177
178 _p1drf = NULL;
179 _p1drb = NULL;
180 _pndrf = NULL;
181 _pndrb = NULL;
182}
183
184
185/* --Methode-- */
186FFTWServer::~FFTWServer()
187{
188 if (_p1df) delete _p1df ;
189 if (_p1db) delete _p1db ;
190 if (_pndf) delete _pndf ;
191 if (_pndb) delete _pndb ;
192
193 if (_p1drf) delete _p1drf ;
194 if (_p1drb) delete _p1drb ;
195 if (_pndrf) delete _pndrf ;
196 if (_pndrb) delete _pndrb ;
197}
198
199/* --Methode-- */
200FFTServerInterface * FFTWServer::Clone()
201{
202 return (new FFTWServer) ;
203}
204
205/* --Methode-- */
206void
207FFTWServer::FFTForward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out)
208{
209 int rank = ckR8.CheckResize(in, out);
210 if (rank == 1) { // One dimensional transform
211 if (_p1df) _p1df->Recreate(in.Size());
212 else _p1df = new FFTWServerPlan(in.Size(), FFTW_FORWARD, false);
213 fftw_one(_p1df->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
214 }
215 else { // Multi dimensional
216 if (in.NbDimensions() > MAXND_FFTW)
217 throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
218 int sz[MAXND_FFTW];
219 FFTWServerPlan::FillSizes(in, sz);
220 // int k1 = 0;
221 // int k2 = 0;
222 // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
223 // sz[k2] = in.Size(k1); k2++;
224 // }
225 if (_pndf) _pndf->Recreate(in.NbDimensions(), sz);
226 else _pndf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_FORWARD, false);
227 fftwnd_one(_pndf->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
228 }
229 if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
230 return;
231}
232
233/* --Methode-- */
234void FFTWServer::FFTBackward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out)
235{
236 int rank = ckR8.CheckResize(in, out);
237 if (rank == 1) { // One dimensional transform
238 if (_p1db) _p1db->Recreate(in.Size());
239 else _p1db = new FFTWServerPlan(in.Size(), FFTW_BACKWARD, false);
240 fftw_one(_p1db->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
241 }
242 else { // Multi dimensional
243 if (in.NbDimensions() > MAXND_FFTW)
244 throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
245 int sz[MAXND_FFTW];
246 FFTWServerPlan::FillSizes(in, sz);
247 // int k1 = 0;
248 // int k2 = 0;
249 // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
250 // sz[k2] = in.Size(k1); k2++;
251 // }
252 if (_pndb) _pndb->Recreate(in.NbDimensions(), sz);
253 else _pndb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_BACKWARD, false);
254 fftwnd_one(_pndb->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
255 }
256
257 return;
258}
259
260
261void FFTWServer::FFTForward(TArray< r_8 > & in, TArray< complex<r_8> > & out)
262{
263 int rank = ckR8.CheckResize(in, out);
264 if (rank == 1) { // One dimensional transform
265 if (_p1drf) _p1drf->Recreate(in.Size());
266 else _p1drf = new FFTWServerPlan(in.Size(), FFTW_REAL_TO_COMPLEX, true);
267 TArray<r_8> outtemp;
268 outtemp.ReSize(in);
269 rfftw_one(_p1drf->rp, (fftw_real *)(in.Data()) , (fftw_real *)(outtemp.Data()));
270 ReShapetoCompl(outtemp, out);
271 }
272 else { // Multi dimensional
273 if (in.NbDimensions() > MAXND_FFTW)
274 throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) rank > MAXND_FFTW !");
275 int sz[MAXND_FFTW];
276 FFTWServerPlan::FillSizes(in, sz);
277 // int k1 = 0;
278 // int k2 = 0;
279 // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
280 // sz[k2] = in.Size(k1); k2++;
281 // }
282 if (_pndrf) _pndrf->Recreate(in.NbDimensions(), sz);
283 else _pndrf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_REAL_TO_COMPLEX, true);
284 rfftwnd_one_real_to_complex(_pndrf->rpnd, (fftw_real *)(in.Data()) ,
285 (fftw_complex *)(out.Data()) );
286 }
287 if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
288 return;
289
290}
291
292
293
294void FFTWServer::FFTBackward(TArray< complex<r_8> > & in, TArray< r_8 > & out,
295 bool usoutsz)
296{
297
298 int rank = ckR8.CheckResize(in, out, usoutsz);
299 if (rank == 1) { // One dimensional transform
300 TArray<r_8> intemp;
301 intemp.ReSize(out);
302 if (_p1drb) _p1drb->Recreate(out.Size());
303 else _p1drb = new FFTWServerPlan(out.Size(), FFTW_COMPLEX_TO_REAL, true);
304
305 ReShapetoReal(in, intemp, usoutsz);
306 // cerr << " DEBUG-FFTWServer::FFTBackward() in = \n" << in << endl;
307 // cerr << " DEBUG-FFTWServer::FFTBackward() intemp = \n" << intemp << endl;
308 rfftw_one(_p1drb->rp, (fftw_real *)(intemp.Data()) , (fftw_real *)(out.Data()));
309
310 }
311 else { // Multi dimensional
312 if (in.NbDimensions() > MAXND_FFTW)
313 throw ParmError("FFTWServer::FFTBackward( complex<r_8>, r_8 ) rank > MAXND_FFTW !");
314 int sz[MAXND_FFTW];
315 FFTWServerPlan::FillSizes(out, sz);
316 // int k1 = 0;
317 // int k2 = 0;
318 // for(k1=out.NbDimensions()-1; k1>=0; k1--) {
319 // sz[k2] = out.Size(k1); k2++;
320 // }
321 if (_pndrb) _pndrb->Recreate(in.NbDimensions(), sz);
322 else _pndrb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_COMPLEX_TO_REAL, true);
323
324 rfftwnd_one_complex_to_real(_pndrb->rpnd, (fftw_complex *)(in.Data()) , (fftw_real *)(out.Data()) );
325 }
326 return;
327}
328
329
330
331/* --Methode-- */
332void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa,
333 bool usz)
334{
335 TVector< complex<r_8> > in(ina);
336 TVector< r_8> out(outa);
337 int n = in.NElts();
338 r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
339 sa_size_t ncs;
340 if (usz) {
341 if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
342 throw SzMismatchError("FFTWServer::ReShapetoReal(..., true) - Wrong output array size ");
343 ncs = out.NElts();
344 }
345 else {
346 ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
347 ncs = 2*n-1 : ncs = 2*n-2;
348
349 if (out.NElts() != ncs) {
350 cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
351 << " out.NElts()= " << out.NElts() << endl;
352 throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
353 }
354 }
355 // if (ncs == 2*n-2) {
356 out(0) = in(0).real();
357 for(int k=1; k<n; k++) {
358 out(k) = in(k).real();
359 out(ncs-k) = in(k).imag();
360 }
361 if (ncs == 2*n-2)
362 out(n-1) = in(n-1).real();
363
364 // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
365}
366
367
368/* --Methode-- */
369void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
370{
371 TVector< r_8> in(ina);
372 TVector< complex<r_8> > out(outa);
373
374 sa_size_t n = in.NElts();
375 sa_size_t ncs = n/2+1;
376 sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
377 if (out.NElts() != ncs)
378 throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
379
380 out(0) = complex<r_8> (in(0),0.);
381 for(int k=1; k<ncs; k++)
382 out(k) = complex<r_8>(in(k), in(n-k));
383 if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
384
385}
386
Note: See TracBrowser for help on using the repository browser.