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

Last change on this file since 1405 was 1405, checked in by ansari, 25 years ago

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