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

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

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