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