Subversion Repositories shark

Rev

Details | 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
 
20
/*
21
 * planner.c -- find the optimal plan
22
 */
23
 
24
/* $Id: rplanner.c,v 1.1.1.1 2002-03-29 14:12:59 pj Exp $ */
25
#ifdef FFTW_USING_CILK
26
#include <cilk.h>
27
#include <cilk-compat.h>
28
#endif
29
 
30
#include <stdlib.h>
31
//#include <stdio.h>
32
 
33
#include <ports/fftw-int.h>
34
#include <ports/rfftw.h>
35
 
36
extern fftw_codelet_desc *rfftw_config[];       /* global from rconfig.c */
37
extern fftw_rgeneric_codelet fftw_hc2hc_forward_generic;
38
extern fftw_rgeneric_codelet fftw_hc2hc_backward_generic;
39
 
40
/* undocumented debugging hook */
41
void (*rfftw_plan_hook)(fftw_plan plan) = (void (*)(fftw_plan)) 0;
42
 
43
/* timing rfftw plans: */
44
static double rfftw_measure_runtime(fftw_plan plan,
45
                                    fftw_real *in, int istride,
46
                                    fftw_real *out, int ostride)
47
{
48
     fftw_time begin, end, start;
49
     double t, tmin;
50
     int i, iter;
51
     int n;
52
     int repeat;
53
 
54
     n = plan->n;
55
 
56
     iter = 1;
57
 
58
     for (;;) {
59
          tmin = 1.0E10;
60
          for (i = 0; i < n; ++i)
61
               in[istride * i] = 0.0;
62
 
63
          start = fftw_get_time();
64
          /* repeat the measurement FFTW_TIME_REPEAT times */
65
          for (repeat = 0; repeat < FFTW_TIME_REPEAT; ++repeat) {
66
               begin = fftw_get_time();
67
               for (i = 0; i < iter; ++i)
68
                    rfftw(plan, 1, in, istride, 0, out, ostride, 0);
69
               end = fftw_get_time();
70
 
71
               t = fftw_time_to_sec(fftw_time_diff(end, begin));
72
               if (t < tmin)
73
                    tmin = t;
74
 
75
               /* do not run for too long */
76
               t = fftw_time_to_sec(fftw_time_diff(end, start));
77
               if (t > FFTW_TIME_LIMIT)
78
                    break;
79
          }
80
 
81
          if (tmin >= FFTW_TIME_MIN)
82
               break;
83
 
84
          iter *= 2;
85
     }
86
 
87
     tmin /= (double) iter;
88
     return tmin;
89
}
90
 
91
/* auxiliary functions */
92
static void rcompute_cost(fftw_plan plan,
93
                          fftw_real *in, int istride,
94
                          fftw_real *out, int ostride)
95
{
96
     if (plan->flags & FFTW_MEASURE)
97
          plan->cost = rfftw_measure_runtime(plan, in, istride, out, ostride);
98
     else {
99
          double c;
100
          c = plan->n * fftw_estimate_node(plan->root);
101
          plan->cost = c;
102
     }
103
}
104
 
105
static void run_plan_hooks(fftw_plan p)
106
{
107
     if (rfftw_plan_hook && p) {
108
          fftw_complete_twiddle(p->root, p->n);
109
          rfftw_plan_hook(p);
110
     }
111
}
112
 
113
/* macrology */
114
#define FOR_ALL_RCODELETS(p) \
115
   fftw_codelet_desc **__q, *p;                         \
116
   for (__q = &rfftw_config[0]; (p = (*__q)); ++__q)
117
 
118
/******************************************
119
 *      Recursive planner                 *
120
 ******************************************/
121
static fftw_plan rplanner(fftw_plan *table, int n,
122
                          fftw_direction dir, int flags,
123
                          fftw_real *, int, fftw_real *, int);
124
 
125
/*
126
 * the planner consists of two parts: one that tries to
127
 * use accumulated wisdom, and one that does not.
128
 * A small driver invokes both parts in sequence
129
 */
130
 
131
/* planner with wisdom: look up the codelet suggested by the wisdom */
132
static fftw_plan rplanner_wisdom(fftw_plan *table, int n,
133
                                 fftw_direction dir, int flags,
134
                                 fftw_real *in, int istride,
135
                                 fftw_real *out, int ostride)
136
{
137
     fftw_plan best = (fftw_plan) 0;
138
     fftw_plan_node *node;
139
     int have_wisdom;
140
     enum fftw_node_type wisdom_type;
141
     int wisdom_signature;
142
 
143
     /* see if we remember any wisdom for this case */
144
     have_wisdom = fftw_wisdom_lookup(n, flags, dir, RFFTW_WISDOM,
145
                                      istride, ostride,
146
                                      &wisdom_type, &wisdom_signature, 0);
147
 
148
     if (!have_wisdom)
149
          return best;
150
 
151
     if (wisdom_type == FFTW_REAL2HC || wisdom_type == FFTW_HC2REAL) {
152
          FOR_ALL_RCODELETS(p) {
153
               if (p->dir == dir && p->type == wisdom_type) {
154
                    /* see if wisdom applies */
155
                    if (wisdom_signature == p->signature &&
156
                        p->size == n) {
157
                         if (wisdom_type == FFTW_REAL2HC)
158
                              node = fftw_make_node_real2hc(n, p);
159
                         else
160
                              node = fftw_make_node_hc2real(n, p);
161
                         best = fftw_make_plan(n, dir, node, flags,
162
                                               p->type, p->signature);
163
                         fftw_use_plan(best);
164
                         run_plan_hooks(best);
165
                         return best;
166
                    }
167
               }
168
          }
169
     }
170
     if (wisdom_type == FFTW_HC2HC) {
171
          FOR_ALL_RCODELETS(p) {
172
               if (p->dir == dir && p->type == wisdom_type) {
173
 
174
                    /* see if wisdom applies */
175
                    if (wisdom_signature == p->signature &&
176
                        p->size > 1 &&
177
                        (n % p->size) == 0) {
178
                         fftw_plan r = rplanner(table, n / p->size, dir,
179
                                                flags, in, istride, out,
180
                                                ostride);
181
                         if (!r)
182
                              continue;
183
                         node = fftw_make_node_hc2hc(n, dir, p,
184
                                                     r->root, flags);
185
                         best = fftw_make_plan(n, dir, node, flags,
186
                                               p->type, p->signature);
187
                         fftw_use_plan(best);
188
                         run_plan_hooks(best);
189
                         fftw_destroy_plan_internal(r);
190
                         return best;
191
                    }
192
               }
193
          }
194
     }
195
     /*
196
      * BUG (or: TODO)  Can we have generic wisdom? This is probably
197
      * an academic question
198
      */
199
 
200
     return best;
201
}
202
 
203
/*
204
 * planner with no wisdom: try all combinations and pick
205
 * the best
206
 */
207
 
208
static fftw_plan rplanner_normal(fftw_plan *table, int n, fftw_direction dir,
209
                                 int flags, fftw_real *in, int istride,
210
                                 fftw_real *out, int ostride)
211
{
212
     fftw_plan best = (fftw_plan) 0;
213
     fftw_plan newplan;
214
     fftw_plan_node *node;
215
 
216
     /* see if we have any codelet that solves the problem */
217
     {
218
          FOR_ALL_RCODELETS(p) {
219
               if (p->dir == dir &&
220
                   (p->type == FFTW_REAL2HC || p->type == FFTW_HC2REAL)) {
221
                    if (p->size == n) {
222
                         if (p->type == FFTW_REAL2HC)
223
                              node = fftw_make_node_real2hc(n, p);
224
                         else
225
                              node = fftw_make_node_hc2real(n, p);
226
                         newplan = fftw_make_plan(n, dir, node, flags,
227
                                                  p->type, p->signature);
228
                         fftw_use_plan(newplan);
229
                         run_plan_hooks(newplan);
230
                         rcompute_cost(newplan, in, istride, out, ostride);
231
                         best = fftw_pick_better(newplan, best);
232
                    }
233
               }
234
          }
235
     }
236
 
237
     /* Then, try all available twiddle codelets */
238
     {
239
          FOR_ALL_RCODELETS(p) {
240
               if (p->dir == dir && p->type == FFTW_HC2HC) {
241
                    if ((n % p->size) == 0 &&
242
                        p->size > 1 &&
243
                        (!best || n != p->size)) {
244
                         fftw_plan r = rplanner(table, n / p->size, dir, flags,
245
                                              in, istride, out, ostride);
246
                         if (!r)
247
                              continue;
248
                         node = fftw_make_node_hc2hc(n, dir, p,
249
                                                     r->root, flags);
250
                         newplan = fftw_make_plan(n, dir, node, flags,
251
                                                  p->type, p->signature);
252
                         fftw_use_plan(newplan);
253
                         run_plan_hooks(newplan);
254
                         fftw_destroy_plan_internal(r);
255
                         rcompute_cost(newplan, in, istride, out, ostride);
256
                         best = fftw_pick_better(newplan, best);
257
                    }
258
               }
259
          }
260
     }
261
 
262
     /*
263
      * Resort to generic codelets for unknown factors, but only if
264
      * n is odd--the rgeneric codelets can't handle even n's.
265
      */
266
     if (n % 2 != 0) {
267
          fftw_rgeneric_codelet *codelet = (dir == FFTW_FORWARD ?
268
                                            fftw_hc2hc_forward_generic :
269
                                            fftw_hc2hc_backward_generic);
270
          int size, prev_size = 0, remaining_factors = n;
271
          fftw_plan r;
272
 
273
          while (remaining_factors > 1) {
274
               size = fftw_factor(remaining_factors);
275
               remaining_factors /= size;
276
 
277
               /* don't try the same factor more than once */
278
               if (size == prev_size)
279
                    continue;
280
               prev_size = size;
281
 
282
               /* Look for codelets corresponding to this factor. */
283
               {
284
                    FOR_ALL_RCODELETS(p) {
285
                         if (p->dir == dir && p->type == FFTW_HC2HC
286
                             && p->size == size) {
287
                              size = 0;
288
                              break;
289
                         }
290
                    }
291
               }
292
 
293
               /*
294
                * only try a generic/rader codelet if there were no
295
                * twiddle codelets for this factor
296
                */
297
               if (!size)
298
                    continue;
299
 
300
               r = rplanner(table, n / size, dir, flags,
301
                            in, istride, out, ostride);
302
 
303
               node = fftw_make_node_rgeneric(n, size, dir, codelet,
304
                                              r->root, flags);
305
               newplan = fftw_make_plan(n, dir, node, flags, FFTW_RGENERIC, 0);
306
               fftw_use_plan(newplan);
307
               run_plan_hooks(newplan);
308
               fftw_destroy_plan_internal(r);
309
               rcompute_cost(newplan, in, istride, out, ostride);
310
               best = fftw_pick_better(newplan, best);
311
          }
312
     }
313
     return best;
314
}
315
 
316
static fftw_plan rplanner(fftw_plan *table, int n, fftw_direction dir,
317
                          int flags, fftw_real *in, int istride,
318
                          fftw_real *out, int ostride)
319
{
320
     fftw_plan best = (fftw_plan) 0;
321
 
322
     /* see if plan has already been computed */
323
     best = fftw_lookup(table, n, flags);
324
     if (best) {
325
          fftw_use_plan(best);
326
          return best;
327
     }
328
     /* try a wise plan */
329
     best = rplanner_wisdom(table, n, dir, flags,
330
                            in, istride, out, ostride);
331
 
332
     if (!best) {
333
          /* No wisdom.  Plan normally. */
334
          best = rplanner_normal(table, n, dir, flags,
335
                                 in, istride, out, ostride);
336
     }
337
     if (best) {
338
          fftw_insert(table, best, n);
339
 
340
          /* remember the wisdom */
341
          fftw_wisdom_add(n, flags, dir, RFFTW_WISDOM,
342
                          istride, ostride,
343
                          best->wisdom_type,
344
                          best->wisdom_signature);
345
     }
346
     return best;
347
}
348
 
349
fftw_plan rfftw_create_plan_specific(int n, fftw_direction dir, int flags,
350
                                     fftw_real *in, int istride,
351
                                     fftw_real *out, int ostride)
352
{
353
     fftw_plan table;
354
     fftw_plan p1;
355
 
356
     /* validate parameters */
357
     if (n <= 0)
358
          return (fftw_plan) 0;
359
 
360
     if ((dir != FFTW_FORWARD) && (dir != FFTW_BACKWARD))
361
          return (fftw_plan) 0;
362
 
363
     fftw_make_empty_table(&table);
364
     p1 = rplanner(&table, n, dir, flags,
365
                   in, istride, out, ostride);
366
     fftw_destroy_table(&table);
367
 
368
     if (p1)
369
          fftw_complete_twiddle(p1->root, n);
370
     return p1;
371
}
372
 
373
fftw_plan rfftw_create_plan(int n, fftw_direction dir, int flags)
374
{
375
     fftw_real *tmp_in;
376
     fftw_real *tmp_out;
377
     fftw_plan p;
378
 
379
     if (flags & FFTW_MEASURE) {
380
          tmp_in = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
381
          tmp_out = (fftw_real *) fftw_malloc(n * sizeof(fftw_real));
382
          if (!tmp_in || !tmp_out)
383
               return 0;
384
 
385
          p = rfftw_create_plan_specific(n, dir, flags,
386
                                         tmp_in, 1, tmp_out, 1);
387
 
388
          fftw_free(tmp_in);
389
          fftw_free(tmp_out);
390
     } else
391
          p = rfftw_create_plan_specific(n, dir, flags,
392
                                 (fftw_real *) 0, 1, (fftw_real *) 0, 1);
393
 
394
     return p;
395
}
396
 
397
void rfftw_destroy_plan(fftw_plan plan)
398
{
399
     fftw_destroy_plan_internal(plan);
400
}
401
/*
402
void rfftw_fprint_plan(FILE *f, fftw_plan p)
403
{
404
     fftw_fprint_plan(f, p);
405
}
406
 
407
void rfftw_print_plan(fftw_plan p)
408
{
409
     rfftw_fprint_plan(stdout, p);
410
}
411
*/