Subversion Repositories shark

Rev

Rev 3 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 1
/*
2
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
3
 *
4
 * This program is free software; you can redistribute it and/or modify
5
 * it under the terms of the GNU General Public License as published by
6
 * the Free Software Foundation; either version 2 of the License, or
7
 * (at your option) any later version.
8
 *
9
 * This program is distributed in the hope that it will be useful,
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12
 * GNU General Public License for more details.
13
 *
14
 * You should have received a copy of the GNU General Public License
15
 * along with this program; if not, write to the Free Software
16
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17
 *
18
 */
19
 
107 pj 20
/* $Id: rfftwnd.c,v 1.2 2003-03-24 11:14:59 pj Exp $ */
2 pj 21
 
107 pj 22
#include <fftw-int.h>
23
#include <rfftw.h>
2 pj 24
 
25
/********************** prototypes for rexec2 routines **********************/
26
 
27
extern void rfftw_real2c_aux(fftw_plan plan, int howmany,
28
                             fftw_real *in, int istride, int idist,
29
                             fftw_complex *out, int ostride, int odist,
30
                             fftw_real *work);
31
extern void rfftw_c2real_aux(fftw_plan plan, int howmany,
32
                             fftw_complex *in, int istride, int idist,
33
                             fftw_real *out, int ostride, int odist,
34
                             fftw_real *work);
35
extern void rfftw_real2c_overlap_aux(fftw_plan plan, int howmany,
36
                                   fftw_real *in, int istride, int idist,
37
                               fftw_complex *out, int ostride, int odist,
38
                                     fftw_real *work);
39
extern void rfftw_c2real_overlap_aux(fftw_plan plan, int howmany,
40
                                fftw_complex *in, int istride, int idist,
41
                                  fftw_real *out, int ostride, int odist,
42
                                     fftw_real *work);
43
 
44
/********************** Initializing the RFFTWND Plan ***********************/
45
 
46
/*
47
 * Create an fftwnd_plan specialized for specific arrays.  (These
48
 * arrays are ignored, however, if they are NULL or if the flags
49
 * do not include FFTW_MEASURE.)  The main advantage of being
50
 * provided arrays like this is that we can do runtime timing
51
 * measurements of our options, without worrying about allocating
52
 * excessive scratch space.
53
 */
54
fftwnd_plan rfftwnd_create_plan_specific(int rank, const int *n,
55
                                         fftw_direction dir, int flags,
56
                                         fftw_real *in, int istride,
57
                                         fftw_real *out, int ostride)
58
{
59
     fftwnd_plan p;
60
     int i;
61
     int rflags = flags & ~FFTW_IN_PLACE;
62
     /* note that we always do rfftw transforms out-of-place in rexec2.c */
63
 
64
     if (flags & FFTW_IN_PLACE) {
65
          out = NULL;
66
          ostride = istride;
67
     }
68
     istride = ostride = 1;     /*
69
                                 * strides don't work yet, since it is not
70
                                 * clear whether they apply to real
71
                                 * or complex data
72
                                 */
73
 
74
     if (!(p = fftwnd_create_plan_aux(rank, n, dir, flags)))
75
          return 0;
76
 
77
     for (i = 0; i < rank - 1; ++i)
78
          p->n_after[i] = (n[rank - 1]/2 + 1) * (p->n_after[i] / n[rank - 1]);
79
     if (rank > 0)
80
          p->n[rank - 1] = n[rank - 1] / 2 + 1;
81
 
82
     p->plans = fftwnd_new_plan_array(rank);
83
     if (rank > 0 && !p->plans) {
84
          rfftwnd_destroy_plan(p);
85
          return 0;
86
     }
87
     if (rank > 0) {
88
          p->plans[rank - 1] = rfftw_create_plan(n[rank - 1], dir, rflags);
89
          if (!p->plans[rank - 1]) {
90
               rfftwnd_destroy_plan(p);
91
               return 0;
92
          }
93
     }
94
     if (rank > 1) {
95
          if (!(flags & FFTW_MEASURE) || in == 0
96
              || (!p->is_in_place && out == 0)) {
97
               if (!fftwnd_create_plans_generic(p->plans, rank - 1, n,
98
                                           dir, flags | FFTW_IN_PLACE)) {
99
                    rfftwnd_destroy_plan(p);
100
                    return 0;
101
               }
102
          } else if (dir == FFTW_COMPLEX_TO_REAL || (flags & FFTW_IN_PLACE)) {
103
               if (!fftwnd_create_plans_specific(p->plans, rank - 1, n,
104
                                                 p->n_after,
105
                                              dir, flags | FFTW_IN_PLACE,
106
                                                 (fftw_complex *) in,
107
                                                 istride,
108
                                                 0, 0)) {
109
                    rfftwnd_destroy_plan(p);
110
                    return 0;
111
               }
112
          } else {
113
               if (!fftwnd_create_plans_specific(p->plans, rank - 1, n,
114
                                                 p->n_after,
115
                                              dir, flags | FFTW_IN_PLACE,
116
                                                 (fftw_complex *) out,
117
                                                 ostride,
118
                                                 0, 0)) {
119
                    rfftwnd_destroy_plan(p);
120
                    return 0;
121
               }
122
          }
123
     }
124
     p->nbuffers = 0;
125
     p->nwork = fftwnd_work_size(rank, p->n, flags | FFTW_IN_PLACE,
126
                                 p->nbuffers + 1);
127
     if (p->nwork && !(flags & FFTW_THREADSAFE)) {
128
          p->work = (fftw_complex *) fftw_malloc(p->nwork
129
                                                 * sizeof(fftw_complex));
130
          if (!p->work) {
131
               rfftwnd_destroy_plan(p);
132
               return 0;
133
          }
134
     }
135
     return p;
136
}
137
 
138
fftwnd_plan rfftw2d_create_plan_specific(int nx, int ny,
139
                                         fftw_direction dir, int flags,
140
                                         fftw_real *in, int istride,
141
                                         fftw_real *out, int ostride)
142
{
143
     int n[2];
144
 
145
     n[0] = nx;
146
     n[1] = ny;
147
 
148
     return rfftwnd_create_plan_specific(2, n, dir, flags,
149
                                         in, istride, out, ostride);
150
}
151
 
152
fftwnd_plan rfftw3d_create_plan_specific(int nx, int ny, int nz,
153
                                         fftw_direction dir, int flags,
154
                                         fftw_real *in, int istride,
155
                                         fftw_real *out, int ostride)
156
{
157
     int n[3];
158
 
159
     n[0] = nx;
160
     n[1] = ny;
161
     n[2] = nz;
162
 
163
     return rfftwnd_create_plan_specific(3, n, dir, flags,
164
                                         in, istride, out, ostride);
165
}
166
 
167
/* Create a generic fftwnd plan: */
168
 
169
fftwnd_plan rfftwnd_create_plan(int rank, const int *n,
170
                                fftw_direction dir, int flags)
171
{
172
     return rfftwnd_create_plan_specific(rank, n, dir, flags, 0, 1, 0, 1);
173
}
174
 
175
fftwnd_plan rfftw2d_create_plan(int nx, int ny,
176
                                fftw_direction dir, int flags)
177
{
178
     return rfftw2d_create_plan_specific(nx, ny, dir, flags, 0, 1, 0, 1);
179
}
180
 
181
fftwnd_plan rfftw3d_create_plan(int nx, int ny, int nz,
182
                                fftw_direction dir, int flags)
183
{
184
     return rfftw3d_create_plan_specific(nx, ny, nz, dir, flags, 0, 1, 0, 1);
185
}
186
 
187
/************************ Freeing the RFFTWND Plan ************************/
188
 
189
void rfftwnd_destroy_plan(fftwnd_plan plan)
190
{
191
     fftwnd_destroy_plan(plan);
192
}
193
 
194
/************************ Printing the RFFTWND Plan ************************/
195
/*
196
void rfftwnd_fprint_plan(FILE *f, fftwnd_plan plan)
197
{
198
     fftwnd_fprint_plan(f, plan);
199
}
200
 
201
void rfftwnd_print_plan(fftwnd_plan plan)
202
{
203
     rfftwnd_fprint_plan(stdout, plan);
204
}
205
*/
206
/*********** Computing the N-Dimensional FFT: Auxiliary Routines ************/
207
 
208
void rfftwnd_real2c_aux(fftwnd_plan p, int cur_dim,
209
                        fftw_real *in, int istride,
210
                        fftw_complex *out, int ostride,
211
                        fftw_real *work)
212
{
213
     int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
214
 
215
     if (cur_dim == p->rank - 2) {
216
          /* just do the last dimension directly: */
217
          if (p->is_in_place)
218
               rfftw_real2c_aux(p->plans[p->rank - 1], n,
219
                                in, istride, (n_after * istride) * 2,
220
                                out, istride, n_after * istride,
221
                                work);
222
          else
223
               rfftw_real2c_aux(p->plans[p->rank - 1], n,
224
                         in, istride, p->plans[p->rank - 1]->n * istride,
225
                                out, ostride, n_after * ostride,
226
                                work);
227
     } else {                   /* we have at least two dimensions to go */
228
          int nr = p->plans[p->rank - 1]->n;
229
          int n_after_r = p->is_in_place ? n_after * 2
230
               : nr * (n_after / (nr/2 + 1));
231
          int i;
232
 
233
          /*
234
           * process the subsequent dimensions recursively, in hyperslabs,
235
           * to get maximum locality:
236
           */
237
          for (i = 0; i < n; ++i)
238
               rfftwnd_real2c_aux(p, cur_dim + 1,
239
                                  in + i * n_after_r * istride, istride,
240
                             out + i * n_after * ostride, ostride, work);
241
     }
242
 
243
     /* do the current dimension (in-place): */
244
     fftw(p->plans[cur_dim], n_after,
245
          out, n_after * ostride, ostride,
246
          (fftw_complex *) work, 1, 0);
247
     /* I hate this cast */
248
}
249
 
250
void rfftwnd_c2real_aux(fftwnd_plan p, int cur_dim,
251
                        fftw_complex *in, int istride,
252
                        fftw_real *out, int ostride,
253
                        fftw_real *work)
254
{
255
     int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
256
 
257
     /* do the current dimension (in-place): */
258
     fftw(p->plans[cur_dim], n_after,
259
          in, n_after * istride, istride,
260
          (fftw_complex *) work, 1, 0);
261
 
262
     if (cur_dim == p->rank - 2) {
263
          /* just do the last dimension directly: */
264
          if (p->is_in_place)
265
               rfftw_c2real_aux(p->plans[p->rank - 1], n,
266
                                in, istride, n_after * istride,
267
                                out, istride, (n_after * istride) * 2,
268
                                work);
269
          else
270
               rfftw_c2real_aux(p->plans[p->rank - 1], n,
271
                                in, istride, n_after * istride,
272
                        out, ostride, p->plans[p->rank - 1]->n * ostride,
273
                                work);
274
     } else {                   /* we have at least two dimensions to go */
275
          int nr = p->plans[p->rank - 1]->n;
276
          int n_after_r = p->is_in_place ? n_after * 2 :
277
               nr * (n_after / (nr/2 + 1));
278
          int i;
279
 
280
          /*
281
           * process the subsequent dimensions recursively, in hyperslabs,
282
           * to get maximum locality:
283
           */
284
          for (i = 0; i < n; ++i)
285
               rfftwnd_c2real_aux(p, cur_dim + 1,
286
                                  in + i * n_after * istride, istride,
287
                           out + i * n_after_r * ostride, ostride, work);
288
     }
289
}
290
 
291
/*
292
 * alternate version of rfftwnd_aux -- this version pushes the howmany
293
 * loop down to the leaves of the computation, for greater locality
294
 * in cases where dist < stride.  It is also required for correctness
295
 * if in==out, and we must call a special version of the executor.
296
 * Note that work must point to 'howmany' copies of its data
297
 * if in == out.
298
 */
299
 
300
void rfftwnd_real2c_aux_howmany(fftwnd_plan p, int cur_dim,
301
                                int howmany,
302
                                fftw_real *in, int istride, int idist,
303
                                fftw_complex *out, int ostride, int odist,
304
                                fftw_complex *work)
305
{
306
     int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
307
     int k;
308
 
309
     if (cur_dim == p->rank - 2) {
310
          /* just do the last dimension directly: */
311
          if (p->is_in_place)
312
               for (k = 0; k < n; ++k)
313
                    rfftw_real2c_overlap_aux(p->plans[p->rank - 1], howmany,
314
                                        in + (k * n_after * istride) * 2,
315
                                             istride, idist,
316
                                           out + (k * n_after * ostride),
317
                                             ostride, odist,
318
                                             (fftw_real *) work);
319
          else {
320
               int nlast = p->plans[p->rank - 1]->n;
321
               for (k = 0; k < n; ++k)
322
                    rfftw_real2c_aux(p->plans[p->rank - 1], howmany,
323
                                     in + k * nlast * istride,
324
                                     istride, idist,
325
                                     out + k * n_after * ostride,
326
                                     ostride, odist,
327
                                     (fftw_real *) work);
328
          }
329
     } else {                   /* we have at least two dimensions to go */
330
          int nr = p->plans[p->rank - 1]->n;
331
          int n_after_r = p->is_in_place ? n_after * 2 :
332
               nr * (n_after / (nr/2 + 1));
333
          int i;
334
 
335
          /*
336
           * process the subsequent dimensions recursively, in hyperslabs,
337
           * to get maximum locality:
338
           */
339
          for (i = 0; i < n; ++i)
340
               rfftwnd_real2c_aux_howmany(p, cur_dim + 1, howmany,
341
                            in + i * n_after_r * istride, istride, idist,
342
                             out + i * n_after * ostride, ostride, odist,
343
                                          work);
344
     }
345
 
346
     /* do the current dimension (in-place): */
347
     for (k = 0; k < n_after; ++k)
348
          fftw(p->plans[cur_dim], howmany,
349
               out + k * ostride, n_after * ostride, odist,
350
               work, 1, 0);
351
}
352
 
353
void rfftwnd_c2real_aux_howmany(fftwnd_plan p, int cur_dim,
354
                                int howmany,
355
                                fftw_complex *in, int istride, int idist,
356
                                fftw_real *out, int ostride, int odist,
357
                                fftw_complex *work)
358
{
359
     int n_after = p->n_after[cur_dim], n = p->n[cur_dim];
360
     int k;
361
 
362
     /* do the current dimension (in-place): */
363
     for (k = 0; k < n_after; ++k)
364
          fftw(p->plans[cur_dim], howmany,
365
               in + k * istride, n_after * istride, idist,
366
               work, 1, 0);
367
 
368
     if (cur_dim == p->rank - 2) {
369
          /* just do the last dimension directly: */
370
          if (p->is_in_place)
371
               for (k = 0; k < n; ++k)
372
                    rfftw_c2real_overlap_aux(p->plans[p->rank - 1], howmany,
373
                                             in + (k * n_after * istride),
374
                                             istride, idist,
375
                                       out + (k * n_after * ostride) * 2,
376
                                             ostride, odist,
377
                                             (fftw_real *) work);
378
          else {
379
               int nlast = p->plans[p->rank - 1]->n;
380
               for (k = 0; k < n; ++k)
381
                    rfftw_c2real_aux(p->plans[p->rank - 1], howmany,
382
                                     in + k * n_after * istride,
383
                                     istride, idist,
384
                                     out + k * nlast * ostride,
385
                                     ostride, odist,
386
                                     (fftw_real *) work);
387
          }
388
     } else {                   /* we have at least two dimensions to go */
389
          int nr = p->plans[p->rank - 1]->n;
390
          int n_after_r = p->is_in_place ? n_after * 2
391
               : nr * (n_after / (nr/2 + 1));
392
          int i;
393
 
394
          /*
395
           * process the subsequent dimensions recursively, in hyperslabs,
396
           * to get maximum locality:
397
           */
398
          for (i = 0; i < n; ++i)
399
               rfftwnd_c2real_aux_howmany(p, cur_dim + 1, howmany,
400
                              in + i * n_after * istride, istride, idist,
401
                           out + i * n_after_r * ostride, ostride, odist,
402
                                          work);
403
     }
404
}
405
 
406
/********** Computing the N-Dimensional FFT: User-Visible Routines **********/
407
 
408
void rfftwnd_real_to_complex(fftwnd_plan p, int howmany,
409
                             fftw_real *in, int istride, int idist,
410
                             fftw_complex *out, int ostride, int odist)
411
{
412
     fftw_complex *work = p->work;
413
     int rank = p->rank;
414
     int free_work = 0;
415
 
416
     if (p->dir != FFTW_REAL_TO_COMPLEX)
417
          fftw_die("rfftwnd_real_to_complex with complex-to-real plan");
418
 
419
#ifdef FFTW_DEBUG
420
     if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE)
421
         && p->nwork && p->work)
422
          fftw_die("bug with FFTW_THREADSAFE flag");
423
#endif
424
 
425
     if (p->is_in_place) {
426
          ostride = istride;
427
          odist = (idist == 1) ? 1 : (idist / 2);       /* ugh */
428
          out = (fftw_complex *) in;
429
          if (howmany > 1 && istride > idist && rank > 0) {
430
               int new_nwork;
431
 
432
               new_nwork = p->n[rank - 1] * howmany;
433
               if (new_nwork > p->nwork) {
434
                    work = (fftw_complex *)
435
                        fftw_malloc(sizeof(fftw_complex) * new_nwork);
436
                    if (!work)
437
                         fftw_die("error allocating work array");
438
                    free_work = 1;
439
               }
440
          }
441
     }
442
     if (p->nwork && !work) {
443
          work = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * p->nwork);
444
          free_work = 1;
445
     }
446
     switch (rank) {
447
         case 0:
448
              break;
449
         case 1:
450
              if (p->is_in_place && howmany > 1 && istride > idist)
451
                   rfftw_real2c_overlap_aux(p->plans[0], howmany,
452
                                            in, istride, idist,
453
                                            out, ostride, odist,
454
                                            (fftw_real *) work);
455
              else
456
                   rfftw_real2c_aux(p->plans[0], howmany,
457
                                    in, istride, idist,
458
                                    out, ostride, odist,
459
                                    (fftw_real *) work);
460
              break;
461
         default:               /* rank >= 2 */
462
              {
463
                   if (howmany > 1 && ostride > odist)
464
                        rfftwnd_real2c_aux_howmany(p, 0, howmany,
465
                                                   in, istride, idist,
466
                                                   out, ostride, odist,
467
                                                   work);
468
                   else {
469
                        int i;
470
 
471
                        for (i = 0; i < howmany; ++i)
472
                             rfftwnd_real2c_aux(p, 0,
473
                                                in + i * idist, istride,
474
                                                out + i * odist, ostride,
475
                                                (fftw_real *) work);
476
                   }
477
              }
478
     }
479
 
480
     if (free_work)
481
          fftw_free(work);
482
}
483
 
484
void rfftwnd_complex_to_real(fftwnd_plan p, int howmany,
485
                             fftw_complex *in, int istride, int idist,
486
                             fftw_real *out, int ostride, int odist)
487
{
488
     fftw_complex *work = p->work;
489
     int rank = p->rank;
490
     int free_work = 0;
491
 
492
     if (p->dir != FFTW_COMPLEX_TO_REAL)
493
          fftw_die("rfftwnd_complex_to_real with real-to-complex plan");
494
 
495
#ifdef FFTW_DEBUG
496
     if (p->rank > 0 && (p->plans[0]->flags & FFTW_THREADSAFE)
497
         && p->nwork && p->work)
498
          fftw_die("bug with FFTW_THREADSAFE flag");
499
#endif
500
 
501
     if (p->is_in_place) {
502
          ostride = istride;
503
          odist = idist;
504
          odist = (idist == 1) ? 1 : (idist * 2);       /* ugh */
505
          out = (fftw_real *) in;
506
          if (howmany > 1 && istride > idist && rank > 0) {
507
               int new_nwork = p->n[rank - 1] * howmany;
508
               if (new_nwork > p->nwork) {
509
                    work = (fftw_complex *)
510
                        fftw_malloc(sizeof(fftw_complex) * new_nwork);
511
                    if (!work)
512
                         fftw_die("error allocating work array");
513
                    free_work = 1;
514
               }
515
          }
516
     }
517
     if (p->nwork && !work) {
518
          work = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * p->nwork);
519
          free_work = 1;
520
     }
521
     switch (rank) {
522
         case 0:
523
              break;
524
         case 1:
525
              if (p->is_in_place && howmany > 1 && istride > idist)
526
                   rfftw_c2real_overlap_aux(p->plans[0], howmany,
527
                                            in, istride, idist,
528
                                            out, ostride, odist,
529
                                            (fftw_real *) work);
530
              else
531
                   rfftw_c2real_aux(p->plans[0], howmany,
532
                                    in, istride, idist,
533
                                    out, ostride, odist,
534
                                    (fftw_real *) work);
535
              break;
536
         default:               /* rank >= 2 */
537
              {
538
                   if (howmany > 1 && ostride > odist)
539
                        rfftwnd_c2real_aux_howmany(p, 0, howmany,
540
                                                   in, istride, idist,
541
                                                   out, ostride, odist,
542
                                                   work);
543
                   else {
544
                        int i;
545
 
546
                        for (i = 0; i < howmany; ++i)
547
                             rfftwnd_c2real_aux(p, 0,
548
                                                in + i * idist, istride,
549
                                                out + i * odist, ostride,
550
                                                (fftw_real *) work);
551
                   }
552
              }
553
     }
554
 
555
     if (free_work)
556
          fftw_free(work);
557
}
558
 
559
void rfftwnd_one_real_to_complex(fftwnd_plan p,
560
                                 fftw_real *in, fftw_complex *out)
561
{
562
     rfftwnd_real_to_complex(p, 1, in, 1, 1, out, 1, 1);
563
}
564
 
565
void rfftwnd_one_complex_to_real(fftwnd_plan p,
566
                                 fftw_complex *in, fftw_real *out)
567
{
568
     rfftwnd_complex_to_real(p, 1, in, 1, 1, out, 1, 1);
569
}