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

Last change on this file since 2791 was 2615, checked in by cmv, 21 years ago

using namespace sophya enleve de machdefs.h, nouveau sopnamsp.h cmv 10/09/2004
File size: 11.6 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);
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{
365 TVector< complex<r_8> > in(ina);
366 TVector< r_8> out(outa);
367 int n = in.NElts();
368 r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
369 sa_size_t ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
370 ncs = 2*n-1 : ncs = 2*n-2;
371
372 if (out.NElts() != ncs) {
373 cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
374 << " out.NElts()= " << out.NElts() << endl;
375 throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
376 }
377 // if (ncs == 2*n-2) {
378 out(0) = in(0).real();
379 for(int k=1; k<n; k++) {
380 out(k) = in(k).real();
381 out(ncs-k) = in(k).imag();
382 }
383 if (ncs == 2*n-2)
384 out(ncs-1) = in(n-1).real();
385
386 // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
387}
388
389
390/* --Methode-- */
391void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
392{
393 TVector< r_8> in(ina);
394 TVector< complex<r_8> > out(outa);
395
396 sa_size_t n = in.NElts();
397 sa_size_t ncs = n/2+1;
398 sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
399 if (out.NElts() != ncs)
400 throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
401
402 out(0) = complex<r_8> (in(0),0.);
403 for(int k=1; k<ncs; k++)
404 out(k) = complex<r_8>(in(k), in(n-k));
405 if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
406
407}
408
Note: See TracBrowser for help on using the repository browser.