source: Sophya/trunk/SophyaExt/IFFTW/fftwserver.cc@ 2994

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

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