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

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

Debug FFTWServer (Pas fini) - Reza 13/2/2001

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