Subversion Repositories shark

Rev

Rev 3 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 1
 
2
/*
3
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
4
 *
5
 * This program is free software; you can redistribute it and/or modify
6
 * it under the terms of the GNU General Public License as published by
7
 * the Free Software Foundation; either version 2 of the License, or
8
 * (at your option) any later version.
9
 *
10
 * This program is distributed in the hope that it will be useful,
11
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13
 * GNU General Public License for more details.
14
 *
15
 * You should have received a copy of the GNU General Public License
16
 * along with this program; if not, write to the Free Software
17
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
18
 *
19
 */
20
 
21
/*
22
 * putils.c -- plan utilities shared by planner.c and rplanner.c
23
 */
24
 
107 pj 25
/* $Id: putils.c,v 1.2 2003-03-24 11:14:54 pj Exp $ */
2 pj 26
#ifdef FFTW_USING_CILK
107 pj 27
#include <cilk.h>
28
#include <cilk-compat.h>
2 pj 29
#endif
30
 
107 pj 31
#include <fftw-int.h>
2 pj 32
#include <stdlib.h>
33
//#include <ports/stdio.h>
34
 
35
int fftw_node_cnt = 0;
36
int fftw_plan_cnt = 0;
37
 
38
/*
39
 * These two constants are used for the FFTW_ESTIMATE flag to help
40
 * create a heuristic plan.  They don't affect FFTW_MEASURE.
41
 */
42
#define NOTW_OPTIMAL_SIZE 32
43
#define TWIDDLE_OPTIMAL_SIZE 12
44
 
45
#define IS_POWER_OF_TWO(n) ((n & (n - 1)) == 0)
46
 
47
/* constructors --- I wish I had ML */
48
fftw_plan_node *fftw_make_node(void)
49
{
50
     fftw_plan_node *p = (fftw_plan_node *)
51
     fftw_malloc(sizeof(fftw_plan_node));
52
     p->refcnt = 0;
53
     fftw_node_cnt++;
54
     return p;
55
}
56
 
57
void fftw_use_node(fftw_plan_node *p)
58
{
59
     ++p->refcnt;
60
}
61
 
62
fftw_plan_node *fftw_make_node_notw(int size, const fftw_codelet_desc *config)
63
{
64
     fftw_plan_node *p = fftw_make_node();
65
 
66
     p->type = config->type;
67
     p->nodeu.notw.size = size;
68
     p->nodeu.notw.codelet = (fftw_notw_codelet *) config->codelet;
69
     p->nodeu.notw.codelet_desc = config;
70
     return p;
71
}
72
 
73
fftw_plan_node *fftw_make_node_real2hc(int size,
74
                                       const fftw_codelet_desc *config)
75
{
76
     fftw_plan_node *p = fftw_make_node();
77
 
78
     p->type = config->type;
79
     p->nodeu.real2hc.size = size;
80
     p->nodeu.real2hc.codelet = (fftw_real2hc_codelet *) config->codelet;
81
     p->nodeu.real2hc.codelet_desc = config;
82
     return p;
83
}
84
 
85
fftw_plan_node *fftw_make_node_hc2real(int size,
86
                                       const fftw_codelet_desc *config)
87
{
88
     fftw_plan_node *p = fftw_make_node();
89
 
90
     p->type = config->type;
91
     p->nodeu.hc2real.size = size;
92
     p->nodeu.hc2real.codelet = (fftw_hc2real_codelet *) config->codelet;
93
     p->nodeu.hc2real.codelet_desc = config;
94
     return p;
95
}
96
 
97
fftw_plan_node *fftw_make_node_twiddle(int n,
98
                                       const fftw_codelet_desc *config,
99
                                       fftw_plan_node *recurse,
100
                                       int flags)
101
{
102
     fftw_plan_node *p = fftw_make_node();
103
 
104
     p->type = config->type;
105
     p->nodeu.twiddle.size = config->size;
106
     p->nodeu.twiddle.codelet = (fftw_twiddle_codelet *) config->codelet;
107
     p->nodeu.twiddle.recurse = recurse;
108
     p->nodeu.twiddle.codelet_desc = config;
109
     fftw_use_node(recurse);
110
     if (flags & FFTW_MEASURE)
111
          p->nodeu.twiddle.tw = fftw_create_twiddle(n, config);
112
     else
113
          p->nodeu.twiddle.tw = 0;
114
     return p;
115
}
116
 
117
fftw_plan_node *fftw_make_node_hc2hc(int n, fftw_direction dir,
118
                                     const fftw_codelet_desc *config,
119
                                     fftw_plan_node *recurse,
120
                                     int flags)
121
{
122
     fftw_plan_node *p = fftw_make_node();
123
 
124
     p->type = config->type;
125
     p->nodeu.hc2hc.size = config->size;
126
     p->nodeu.hc2hc.dir = dir;
127
     p->nodeu.hc2hc.codelet = (fftw_hc2hc_codelet *) config->codelet;
128
     p->nodeu.hc2hc.recurse = recurse;
129
     p->nodeu.hc2hc.codelet_desc = config;
130
     fftw_use_node(recurse);
131
     if (flags & FFTW_MEASURE)
132
          p->nodeu.hc2hc.tw = fftw_create_twiddle(n, config);
133
     else
134
          p->nodeu.hc2hc.tw = 0;
135
     return p;
136
}
137
 
138
fftw_plan_node *fftw_make_node_generic(int n, int size,
139
                                       fftw_generic_codelet *codelet,
140
                                       fftw_plan_node *recurse,
141
                                       int flags)
142
{
143
     fftw_plan_node *p = fftw_make_node();
144
 
145
     p->type = FFTW_GENERIC;
146
     p->nodeu.generic.size = size;
147
     p->nodeu.generic.codelet = codelet;
148
     p->nodeu.generic.recurse = recurse;
149
     fftw_use_node(recurse);
150
 
151
     if (flags & FFTW_MEASURE)
152
          p->nodeu.generic.tw = fftw_create_twiddle(n,
153
                                          (const fftw_codelet_desc *) 0);
154
     else
155
          p->nodeu.generic.tw = 0;
156
     return p;
157
}
158
 
159
fftw_plan_node *fftw_make_node_rgeneric(int n, int size,
160
                                        fftw_direction dir,
161
                                        fftw_rgeneric_codelet *codelet,
162
                                        fftw_plan_node *recurse,
163
                                        int flags)
164
{
165
     fftw_plan_node *p = fftw_make_node();
166
 
167
     if (size % 2 == 0 || (n / size) % 2 == 0)
168
          fftw_die("invalid size for rgeneric codelet\n");
169
 
170
     p->type = FFTW_RGENERIC;
171
     p->nodeu.rgeneric.size = size;
172
     p->nodeu.rgeneric.dir = dir;
173
     p->nodeu.rgeneric.codelet = codelet;
174
     p->nodeu.rgeneric.recurse = recurse;
175
     fftw_use_node(recurse);
176
 
177
     if (flags & FFTW_MEASURE)
178
          p->nodeu.rgeneric.tw = fftw_create_twiddle(n,
179
                                          (const fftw_codelet_desc *) 0);
180
     else
181
          p->nodeu.rgeneric.tw = 0;
182
     return p;
183
}
184
 
185
/*
186
 * Note that these two Rader-related things must go here, rather than
187
 * in rader.c, in order that putils.c (and rplanner.c) won't depend
188
 * upon rader.c.
189
 */
190
 
191
fftw_rader_data *fftw_rader_top = NULL;
192
 
193
static void fftw_destroy_rader(fftw_rader_data * d)
194
{
195
     if (d) {
196
          d->refcount--;
197
          if (d->refcount <= 0) {
198
               fftw_rader_data *cur = fftw_rader_top, *prev = NULL;
199
 
200
               while (cur && cur != d) {
201
                    prev = cur;
202
                    cur = cur->next;
203
               }
204
               if (!cur)
205
                    fftw_die("invalid Rader data pointer\n");
206
 
207
               if (prev)
208
                    prev->next = d->next;
209
               else
210
                    fftw_rader_top = d->next;
211
 
212
               fftw_destroy_plan_internal(d->plan);
213
               fftw_free(d->omega);
214
               fftw_free(d->cdesc);
215
               fftw_free(d);
216
          }
217
     }
218
}
219
 
220
static void destroy_tree(fftw_plan_node *p)
221
{
222
     if (p) {
223
          --p->refcnt;
224
          if (p->refcnt == 0) {
225
               switch (p->type) {
226
                   case FFTW_NOTW:
227
                   case FFTW_REAL2HC:
228
                   case FFTW_HC2REAL:
229
                        break;
230
 
231
                   case FFTW_TWIDDLE:
232
                        if (p->nodeu.twiddle.tw)
233
                             fftw_destroy_twiddle(p->nodeu.twiddle.tw);
234
                        destroy_tree(p->nodeu.twiddle.recurse);
235
                        break;
236
 
237
                   case FFTW_HC2HC:
238
                        if (p->nodeu.hc2hc.tw)
239
                             fftw_destroy_twiddle(p->nodeu.hc2hc.tw);
240
                        destroy_tree(p->nodeu.hc2hc.recurse);
241
                        break;
242
 
243
                   case FFTW_GENERIC:
244
                        if (p->nodeu.generic.tw)
245
                             fftw_destroy_twiddle(p->nodeu.generic.tw);
246
                        destroy_tree(p->nodeu.generic.recurse);
247
                        break;
248
 
249
                   case FFTW_RADER:
250
                        if (p->nodeu.rader.tw)
251
                             fftw_destroy_twiddle(p->nodeu.rader.tw);
252
                        if (p->nodeu.rader.rader_data)
253
                             fftw_destroy_rader(p->nodeu.rader.rader_data);
254
                        destroy_tree(p->nodeu.rader.recurse);
255
                        break;
256
 
257
                   case FFTW_RGENERIC:
258
                        if (p->nodeu.rgeneric.tw)
259
                             fftw_destroy_twiddle(p->nodeu.rgeneric.tw);
260
                        destroy_tree(p->nodeu.rgeneric.recurse);
261
                        break;
262
               }
263
 
264
               fftw_free(p);
265
               fftw_node_cnt--;
266
          }
267
     }
268
}
269
 
270
/* create a plan with twiddle factors, and other bells and whistles */
271
fftw_plan fftw_make_plan(int n, fftw_direction dir,
272
                         fftw_plan_node *root, int flags,
273
                         enum fftw_node_type wisdom_type,
274
                         int wisdom_signature)
275
{
276
     fftw_plan p = (fftw_plan) fftw_malloc(sizeof(struct fftw_plan_struct));
277
 
278
     p->n = n;
279
     p->dir = dir;
280
     p->flags = flags;
281
     fftw_use_node(root);
282
     p->root = root;
283
     p->cost = 0.0;
284
     p->wisdom_type = wisdom_type;
285
     p->wisdom_signature = wisdom_signature;
286
     p->next = (fftw_plan) 0;
287
     p->refcnt = 0;
288
     fftw_plan_cnt++;
289
     return p;
290
}
291
 
292
/*
293
 * complete with twiddle factors (because nodes don't have
294
 * them when FFTW_ESTIMATE is set)
295
 */
296
void fftw_complete_twiddle(fftw_plan_node *p, int n)
297
{
298
     int r;
299
     switch (p->type) {
300
         case FFTW_NOTW:
301
         case FFTW_REAL2HC:
302
         case FFTW_HC2REAL:
303
              break;
304
 
305
         case FFTW_TWIDDLE:
306
              r = p->nodeu.twiddle.size;
307
              if (!p->nodeu.twiddle.tw)
308
                   p->nodeu.twiddle.tw =
309
                       fftw_create_twiddle(n, p->nodeu.twiddle.codelet_desc);
310
              fftw_complete_twiddle(p->nodeu.twiddle.recurse, n / r);
311
              break;
312
 
313
         case FFTW_HC2HC:
314
              r = p->nodeu.hc2hc.size;
315
              if (!p->nodeu.hc2hc.tw)
316
                   p->nodeu.hc2hc.tw =
317
                       fftw_create_twiddle(n, p->nodeu.hc2hc.codelet_desc);
318
              fftw_complete_twiddle(p->nodeu.hc2hc.recurse, n / r);
319
              break;
320
 
321
         case FFTW_GENERIC:
322
              r = p->nodeu.generic.size;
323
              if (!p->nodeu.generic.tw)
324
                   p->nodeu.generic.tw =
325
                       fftw_create_twiddle(n, (const fftw_codelet_desc *) 0);
326
              fftw_complete_twiddle(p->nodeu.generic.recurse, n / r);
327
              break;
328
 
329
         case FFTW_RADER:
330
              r = p->nodeu.rader.size;
331
              if (!p->nodeu.rader.tw)
332
                   p->nodeu.rader.tw =
333
                       fftw_create_twiddle(n, p->nodeu.rader.rader_data->cdesc);
334
              fftw_complete_twiddle(p->nodeu.rader.recurse, n / r);
335
              break;
336
 
337
         case FFTW_RGENERIC:
338
              r = p->nodeu.rgeneric.size;
339
              if (!p->nodeu.rgeneric.tw)
340
                   p->nodeu.rgeneric.tw =
341
                       fftw_create_twiddle(n, (const fftw_codelet_desc *) 0);
342
              fftw_complete_twiddle(p->nodeu.rgeneric.recurse, n / r);
343
              break;
344
 
345
     }
346
}
347
 
348
void fftw_use_plan(fftw_plan p)
349
{
350
     ++p->refcnt;
351
}
352
 
353
void fftw_destroy_plan_internal(fftw_plan p)
354
{
355
     --p->refcnt;
356
 
357
     if (p->refcnt == 0) {
358
          destroy_tree(p->root);
359
          fftw_plan_cnt--;
360
          fftw_free(p);
361
     }
362
}
363
 
364
/* end of constructors */
365
 
366
/* management of plan tables */
367
void fftw_make_empty_table(fftw_plan *table)
368
{
369
     *table = (fftw_plan) 0;
370
}
371
 
372
void fftw_insert(fftw_plan *table, fftw_plan this_plan, int n)
373
{
374
     fftw_use_plan(this_plan);
375
     this_plan->n = n;
376
     this_plan->next = *table;
377
     *table = this_plan;
378
}
379
 
380
fftw_plan fftw_lookup(fftw_plan *table, int n, int flags)
381
{
382
     fftw_plan p;
383
 
384
     for (p = *table; p &&
385
          ((p->n != n) || (p->flags != flags)); p = p->next);
386
 
387
     return p;
388
}
389
 
390
void fftw_destroy_table(fftw_plan *table)
391
{
392
     fftw_plan p, q;
393
 
394
     for (p = *table; p; p = q) {
395
          q = p->next;
396
          fftw_destroy_plan_internal(p);
397
     }
398
}
399
 
400
double fftw_estimate_node(fftw_plan_node *p)
401
{
402
     int k;
403
 
404
     switch (p->type) {
405
         case FFTW_NOTW:
406
              k = p->nodeu.notw.size;
407
              goto common1;
408
 
409
         case FFTW_REAL2HC:
410
              k = p->nodeu.real2hc.size;
411
              goto common1;
412
 
413
         case FFTW_HC2REAL:
414
              k = p->nodeu.hc2real.size;
415
            common1:
416
              return 1.0 + 0.1 * (k - NOTW_OPTIMAL_SIZE) *
417
                  (k - NOTW_OPTIMAL_SIZE);
418
 
419
         case FFTW_TWIDDLE:
420
              k = p->nodeu.twiddle.size;
421
              return 1.0 + 0.1 * (k - TWIDDLE_OPTIMAL_SIZE) *
422
                  (k - TWIDDLE_OPTIMAL_SIZE)
423
                  + fftw_estimate_node(p->nodeu.twiddle.recurse);
424
 
425
         case FFTW_HC2HC:
426
              k = p->nodeu.hc2hc.size;
427
              return 1.0 + 0.1 * (k - TWIDDLE_OPTIMAL_SIZE) *
428
                  (k - TWIDDLE_OPTIMAL_SIZE)
429
                  + fftw_estimate_node(p->nodeu.hc2hc.recurse);
430
 
431
         case FFTW_GENERIC:
432
              k = p->nodeu.generic.size;
433
              return 10.0 + k * k
434
                  + fftw_estimate_node(p->nodeu.generic.recurse);
435
 
436
         case FFTW_RADER:
437
              k = p->nodeu.rader.size;
438
              return 10.0 + 10 * k
439
                  + fftw_estimate_node(p->nodeu.rader.recurse);
440
 
441
         case FFTW_RGENERIC:
442
              k = p->nodeu.rgeneric.size;
443
              return 10.0 + k * k
444
                  + fftw_estimate_node(p->nodeu.rgeneric.recurse);
445
     }
446
     return 1.0E20;
447
}
448
 
449
/* pick the better of two plans and destroy the other one. */
450
fftw_plan fftw_pick_better(fftw_plan p1, fftw_plan p2)
451
{
452
     if (!p1)
453
          return p2;
454
 
455
     if (!p2)
456
          return p1;
457
 
458
     if (p1->cost > p2->cost) {
459
          fftw_destroy_plan_internal(p1);
460
          return p2;
461
     } else {
462
          fftw_destroy_plan_internal(p2);
463
          return p1;
464
     }
465
}
466
 
467
/* find the smallest prime factor of n */
468
int fftw_factor(int n)
469
{
470
     int r;
471
 
472
     /* try 2 */
473
     if ((n & 1) == 0)
474
          return 2;
475
 
476
     /* try odd numbers up to sqrt(n) */
477
     for (r = 3; r * r <= n; r += 2)
478
          if (n % r == 0)
479
               return r;
480
 
481
     /* n is prime */
482
     return n;
483
}
484
/*
485
static void print_node(FILE *f, fftw_plan_node *p, int indent)
486
{
487
     if (p) {
488
          switch (p->type) {
489
              case FFTW_NOTW:
490
                   fprintf(f, "%*sFFTW_NOTW %d\n", indent, "",
491
                           p->nodeu.notw.size);
492
                   break;
493
              case FFTW_REAL2HC:
494
                   fprintf(f, "%*sFFTW_REAL2HC %d\n", indent, "",
495
                           p->nodeu.real2hc.size);
496
                   break;
497
              case FFTW_HC2REAL:
498
                   fprintf(f, "%*sFFTW_HC2REAL %d\n", indent, "",
499
                           p->nodeu.hc2real.size);
500
                   break;
501
              case FFTW_TWIDDLE:
502
                   fprintf(f, "%*sFFTW_TWIDDLE %d\n", indent, "",
503
                           p->nodeu.twiddle.size);
504
                   print_node(f, p->nodeu.twiddle.recurse, indent);
505
                   break;
506
              case FFTW_HC2HC:
507
                   fprintf(f, "%*sFFTW_HC2HC %d\n", indent, "",
508
                           p->nodeu.hc2hc.size);
509
                   print_node(f, p->nodeu.hc2hc.recurse, indent);
510
                   break;
511
              case FFTW_GENERIC:
512
                   fprintf(f, "%*sFFTW_GENERIC %d\n", indent, "",
513
                           p->nodeu.generic.size);
514
                   print_node(f, p->nodeu.generic.recurse, indent);
515
                   break;
516
              case FFTW_RADER:
517
                   fprintf(f, "%*sFFTW_RADER %d\n", indent, "",
518
                           p->nodeu.rader.size);
519
 
520
                   fprintf(f, "%*splan for size %d convolution:\n",
521
                           indent + 4, "", p->nodeu.rader.size - 1);
522
                   print_node(f, p->nodeu.rader.rader_data->plan->root,
523
                              indent + 6);
524
 
525
                   print_node(f, p->nodeu.rader.recurse, indent);
526
                   break;
527
              case FFTW_RGENERIC:
528
                   fprintf(f, "%*sFFTW_RGENERIC %d\n", indent, "",
529
                           p->nodeu.rgeneric.size);
530
                   print_node(f, p->nodeu.rgeneric.recurse, indent);
531
                   break;
532
          }
533
     }
534
}
535
 
536
void fftw_fprint_plan(FILE *f, fftw_plan p)
537
{
538
 
539
     fprintf(f, "plan: (cost = %e)\n", p->cost);
540
     print_node(f, p->root, 0);
541
}
542
 
543
void fftw_print_plan(fftw_plan p)
544
{
545
     fftw_fprint_plan(stdout, p);
546
}
547
*/
548
size_t fftw_sizeof_fftw_real(void)
549
{
550
     return(sizeof(fftw_real));
551
}