Branch data Line data Source code
1 : : /* StarPU --- Runtime system for heterogeneous multicore architectures.
2 : : *
3 : : * Copyright (C) 2009-2012 Université de Bordeaux 1
4 : : * Copyright (C) 2010, 2011, 2012, 2013 Centre National de la Recherche Scientifique
5 : : *
6 : : * StarPU is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU Lesser General Public License as published by
8 : : * the Free Software Foundation; either version 2.1 of the License, or (at
9 : : * your option) any later version.
10 : : *
11 : : * StarPU is distributed in the hope that it will be useful, but
12 : : * WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
14 : : *
15 : : * See the GNU Lesser General Public License in COPYING.LGPL for more details.
16 : : */
17 : :
18 : : /*
19 : : *
20 : : * Dumb parallel version
21 : : *
22 : : */
23 : :
24 : : #define DIV_1D 64
25 : :
26 : : /*
27 : : * Overall strategy for an fft of size n:
28 : : * - perform n1 ffts of size n2
29 : : * - twiddle
30 : : * - perform n2 ffts of size n1
31 : : *
32 : : * - n1 defaults to DIV_1D, thus n2 defaults to n / DIV_1D.
33 : : *
34 : : * Precise tasks:
35 : : *
36 : : * - twist1: twist the whole n-element input (called "in") into n1 chunks of
37 : : * size n2, by using n1 tasks taking the whole n-element input as a
38 : : * R parameter and one n2 output as a W parameter. The result is
39 : : * called twisted1.
40 : : * - fft1: perform n1 (n2) ffts, by using n1 tasks doing one fft each. Also
41 : : * twiddle the result to prepare for the fft2. The result is called
42 : : * fft1.
43 : : * - join: depends on all the fft1s, to gather the n1 results of size n2 in
44 : : * the fft1 vector.
45 : : * - twist2: twist the fft1 vector into n2 chunks of size n1, called twisted2.
46 : : * since n2 is typically very large, this step is divided in DIV_1D
47 : : * tasks, each of them performing n2/DIV_1D of them
48 : : * - fft2: perform n2 ffts of size n1. This is divided in DIV_1D tasks of
49 : : * n2/DIV_1D ffts, to be performed in batches. The result is called
50 : : * fft2.
51 : : * - twist3: twist back the result of the fft2s above into the output buffer.
52 : : * Only implemented on CPUs for simplicity of the gathering.
53 : : *
54 : : * The tag space thus uses 3 dimensions:
55 : : * - the number of the plan.
56 : : * - the step (TWIST1, FFT1, JOIN, TWIST2, FFT2, TWIST3, END)
57 : : * - an index i between 0 and DIV_1D-1.
58 : : */
59 : :
60 : : #define STEP_TAG_1D(plan, step, i) _STEP_TAG(plan, step, i)
61 : :
62 : : #ifdef __STARPU_USE_CUDA
63 : : /* twist1:
64 : : *
65 : : * Twist the full input vector (first parameter) into one chunk of size n2
66 : : * (second parameter) */
67 : : static void
68 : 0 : STARPUFFT(twist1_1d_kernel_gpu)(void *descr[], void *_args)
69 : : {
70 : 0 : struct STARPUFFT(args) *args = _args;
71 : 0 : STARPUFFT(plan) plan = args->plan;
72 : 0 : int i = args->i;
73 : 0 : int n1 = plan->n1[0];
74 : 0 : int n2 = plan->n2[0];
75 : :
76 : 0 : _cufftComplex * restrict in = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[0]);
77 : 0 : _cufftComplex * restrict twisted1 = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[1]);
78 : :
79 : 0 : STARPUFFT(cuda_twist1_1d_host)(in, twisted1, i, n1, n2);
80 : :
81 : 0 : cudaStreamSynchronize(starpu_cuda_get_local_stream());
82 : 0 : }
83 : :
84 : : /* fft1:
85 : : *
86 : : * Perform one fft of size n2 */
87 : : static void
88 : 0 : STARPUFFT(fft1_1d_plan_gpu)(void *args)
89 : : {
90 : 0 : STARPUFFT(plan) plan = args;
91 : 0 : int n2 = plan->n2[0];
92 : 0 : int workerid = starpu_worker_get_id();
93 : : cufftResult cures;
94 : :
95 : 0 : cures = cufftPlan1d(&plan->plans[workerid].plan1_cuda, n2, _CUFFT_C2C, 1);
96 [ # # ]: 0 : if (cures != CUFFT_SUCCESS)
97 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
98 : 0 : cufftSetStream(plan->plans[workerid].plan1_cuda, starpu_cuda_get_local_stream());
99 [ # # ]: 0 : if (cures != CUFFT_SUCCESS)
100 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
101 : 0 : }
102 : :
103 : : static void
104 : 0 : STARPUFFT(fft1_1d_kernel_gpu)(void *descr[], void *_args)
105 : : {
106 : 0 : struct STARPUFFT(args) *args = _args;
107 : 0 : STARPUFFT(plan) plan = args->plan;
108 : 0 : int i = args->i;
109 : 0 : int n2 = plan->n2[0];
110 : : cufftResult cures;
111 : :
112 : 0 : _cufftComplex * restrict in = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[0]);
113 : 0 : _cufftComplex * restrict out = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[1]);
114 : 0 : const _cufftComplex * restrict roots = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[2]);
115 : :
116 : 0 : int workerid = starpu_worker_get_id();
117 : :
118 : 0 : task_per_worker[workerid]++;
119 : :
120 [ # # ]: 0 : cures = _cufftExecC2C(plan->plans[workerid].plan1_cuda, in, out, plan->sign == -1 ? CUFFT_FORWARD : CUFFT_INVERSE);
121 [ # # ]: 0 : if (cures != CUFFT_SUCCESS)
122 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
123 : :
124 : 0 : STARPUFFT(cuda_twiddle_1d_host)(out, roots, n2, i);
125 : :
126 : 0 : cudaStreamSynchronize(starpu_cuda_get_local_stream());
127 : 0 : }
128 : :
129 : : /* fft2:
130 : : *
131 : : * Perform n3 = n2/DIV_1D ffts of size n1 */
132 : : static void
133 : 0 : STARPUFFT(fft2_1d_plan_gpu)(void *args)
134 : : {
135 : 0 : STARPUFFT(plan) plan = args;
136 : 0 : int n1 = plan->n1[0];
137 : 0 : int n2 = plan->n2[0];
138 : 0 : int n3 = n2/DIV_1D;
139 : : cufftResult cures;
140 : 0 : int workerid = starpu_worker_get_id();
141 : :
142 : 0 : cures = cufftPlan1d(&plan->plans[workerid].plan2_cuda, n1, _CUFFT_C2C, n3);
143 [ # # ]: 0 : if (cures != CUFFT_SUCCESS)
144 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
145 : 0 : cufftSetStream(plan->plans[workerid].plan2_cuda, starpu_cuda_get_local_stream());
146 [ # # ]: 0 : if (cures != CUFFT_SUCCESS)
147 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
148 : 0 : }
149 : :
150 : : static void
151 : 0 : STARPUFFT(fft2_1d_kernel_gpu)(void *descr[], void *_args)
152 : : {
153 : 0 : struct STARPUFFT(args) *args = _args;
154 : 0 : STARPUFFT(plan) plan = args->plan;
155 : : cufftResult cures;
156 : :
157 : 0 : _cufftComplex * restrict in = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[0]);
158 : 0 : _cufftComplex * restrict out = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[1]);
159 : :
160 : 0 : int workerid = starpu_worker_get_id();
161 : :
162 : 0 : task_per_worker[workerid]++;
163 : :
164 : : /* NOTE using batch support */
165 [ # # ]: 0 : cures = _cufftExecC2C(plan->plans[workerid].plan2_cuda, in, out, plan->sign == -1 ? CUFFT_FORWARD : CUFFT_INVERSE);
166 [ # # ]: 0 : if (cures != CUFFT_SUCCESS)
167 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
168 : :
169 : 0 : cudaStreamSynchronize(starpu_cuda_get_local_stream());
170 : 0 : }
171 : : #endif
172 : :
173 : : /* twist1:
174 : : *
175 : : * Twist the full input vector (first parameter) into one chunk of size n2
176 : : * (second parameter) */
177 : : static void
178 : 0 : STARPUFFT(twist1_1d_kernel_cpu)(void *descr[], void *_args)
179 : : {
180 : 0 : struct STARPUFFT(args) *args = _args;
181 : 0 : STARPUFFT(plan) plan = args->plan;
182 : 0 : int i = args->i;
183 : : int j;
184 : 0 : int n1 = plan->n1[0];
185 : 0 : int n2 = plan->n2[0];
186 : :
187 : 0 : STARPUFFT(complex) * restrict in = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[0]);
188 : 0 : STARPUFFT(complex) * restrict twisted1 = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[1]);
189 : :
190 : : /* printf("twist1 %d %g\n", i, (double) cabs(plan->in[i])); */
191 : :
192 [ # # ]: 0 : for (j = 0; j < n2; j++)
193 : 0 : twisted1[j] = in[i+j*n1];
194 : 0 : }
195 : :
196 : : #ifdef STARPU_HAVE_FFTW
197 : : /* fft1:
198 : : *
199 : : * Perform one fft of size n2 */
200 : : static void
201 : 0 : STARPUFFT(fft1_1d_kernel_cpu)(void *descr[], void *_args)
202 : : {
203 : 0 : struct STARPUFFT(args) *args = _args;
204 : 0 : STARPUFFT(plan) plan = args->plan;
205 : 0 : int i = args->i;
206 : : int j;
207 : 0 : int n2 = plan->n2[0];
208 : 0 : int workerid = starpu_worker_get_id();
209 : :
210 : 0 : task_per_worker[workerid]++;
211 : :
212 : 0 : STARPUFFT(complex) * restrict twisted1 = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[0]);
213 : 0 : STARPUFFT(complex) * restrict fft1 = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[1]);
214 : :
215 : : /* printf("fft1 %d %g\n", i, (double) cabs(twisted1[0])); */
216 : :
217 : 0 : _FFTW(execute_dft)(plan->plans[workerid].plan1_cpu, twisted1, fft1);
218 : :
219 : : /* twiddle fft1 buffer */
220 [ # # ]: 0 : for (j = 0; j < n2; j++)
221 : 0 : fft1[j] = fft1[j] * plan->roots[0][i*j];
222 : 0 : }
223 : : #endif
224 : :
225 : : /* twist2:
226 : : *
227 : : * Twist the full vector (results of the fft1s) into one package of n2/DIV_1D
228 : : * chunks of size n1 */
229 : : static void
230 : 0 : STARPUFFT(twist2_1d_kernel_cpu)(void *descr[], void *_args)
231 : : {
232 : 0 : struct STARPUFFT(args) *args = _args;
233 : 0 : STARPUFFT(plan) plan = args->plan;
234 : 0 : int jj = args->jj; /* between 0 and DIV_1D */
235 : : int jjj; /* beetween 0 and n3 */
236 : : int i;
237 : 0 : int n1 = plan->n1[0];
238 : 0 : int n2 = plan->n2[0];
239 : 0 : int n3 = n2/DIV_1D;
240 : :
241 : 0 : STARPUFFT(complex) * restrict twisted2 = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[0]);
242 : :
243 : : /* printf("twist2 %d %g\n", jj, (double) cabs(plan->fft1[jj])); */
244 : :
245 [ # # ]: 0 : for (jjj = 0; jjj < n3; jjj++) {
246 : 0 : int j = jj * n3 + jjj;
247 [ # # ]: 0 : for (i = 0; i < n1; i++)
248 : 0 : twisted2[jjj*n1+i] = plan->fft1[i*n2+j];
249 : : }
250 : 0 : }
251 : :
252 : : #ifdef STARPU_HAVE_FFTW
253 : : /* fft2:
254 : : *
255 : : * Perform n3 = n2/DIV_1D ffts of size n1 */
256 : : static void
257 : 0 : STARPUFFT(fft2_1d_kernel_cpu)(void *descr[], void *_args)
258 : : {
259 : 0 : struct STARPUFFT(args) *args = _args;
260 : 0 : STARPUFFT(plan) plan = args->plan;
261 : : /* int jj = args->jj; */
262 : 0 : int workerid = starpu_worker_get_id();
263 : :
264 : 0 : task_per_worker[workerid]++;
265 : :
266 : 0 : STARPUFFT(complex) * restrict twisted2 = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[0]);
267 : 0 : STARPUFFT(complex) * restrict fft2 = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[1]);
268 : :
269 : : /* printf("fft2 %d %g\n", jj, (double) cabs(twisted2[plan->totsize4-1])); */
270 : :
271 : 0 : _FFTW(execute_dft)(plan->plans[workerid].plan2_cpu, twisted2, fft2);
272 : 0 : }
273 : : #endif
274 : :
275 : : /* twist3:
276 : : *
277 : : * Spread the package of n2/DIV_1D chunks of size n1 into the output vector */
278 : : static void
279 : 0 : STARPUFFT(twist3_1d_kernel_cpu)(void *descr[], void *_args)
280 : : {
281 : 0 : struct STARPUFFT(args) *args = _args;
282 : 0 : STARPUFFT(plan) plan = args->plan;
283 : 0 : int jj = args->jj; /* between 0 and DIV_1D */
284 : : int jjj; /* beetween 0 and n3 */
285 : : int i;
286 : 0 : int n1 = plan->n1[0];
287 : 0 : int n2 = plan->n2[0];
288 : 0 : int n3 = n2/DIV_1D;
289 : :
290 : 0 : const STARPUFFT(complex) * restrict fft2 = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[0]);
291 : :
292 : : /* printf("twist3 %d %g\n", jj, (double) cabs(fft2[0])); */
293 : :
294 [ # # ]: 0 : for (jjj = 0; jjj < n3; jjj++) {
295 : 0 : int j = jj * n3 + jjj;
296 [ # # ]: 0 : for (i = 0; i < n1; i++)
297 : 0 : plan->out[i*n2+j] = fft2[jjj*n1+i];
298 : : }
299 : 0 : }
300 : :
301 : : /* Performance models for the 5 kinds of tasks */
302 : : static struct starpu_perfmodel STARPUFFT(twist1_1d_model) = {
303 : : .type = STARPU_HISTORY_BASED,
304 : : .symbol = TYPE"twist1_1d"
305 : : };
306 : :
307 : : static struct starpu_perfmodel STARPUFFT(fft1_1d_model) = {
308 : : .type = STARPU_HISTORY_BASED,
309 : : .symbol = TYPE"fft1_1d"
310 : : };
311 : :
312 : : static struct starpu_perfmodel STARPUFFT(twist2_1d_model) = {
313 : : .type = STARPU_HISTORY_BASED,
314 : : .symbol = TYPE"twist2_1d"
315 : : };
316 : :
317 : : static struct starpu_perfmodel STARPUFFT(fft2_1d_model) = {
318 : : .type = STARPU_HISTORY_BASED,
319 : : .symbol = TYPE"fft2_1d"
320 : : };
321 : :
322 : : static struct starpu_perfmodel STARPUFFT(twist3_1d_model) = {
323 : : .type = STARPU_HISTORY_BASED,
324 : : .symbol = TYPE"twist3_1d"
325 : : };
326 : :
327 : : /* codelet pointers for the 5 kinds of tasks */
328 : : static struct starpu_codelet STARPUFFT(twist1_1d_codelet) = {
329 : : .where =
330 : : #ifdef __STARPU_USE_CUDA
331 : : STARPU_CUDA|
332 : : #endif
333 : : STARPU_CPU,
334 : : #ifdef __STARPU_USE_CUDA
335 : : .cuda_funcs = {STARPUFFT(twist1_1d_kernel_gpu), NULL},
336 : : #endif
337 : : .cpu_funcs = {STARPUFFT(twist1_1d_kernel_cpu), NULL},
338 : : CAN_EXECUTE
339 : : .model = &STARPUFFT(twist1_1d_model),
340 : : .nbuffers = 2,
341 : : .modes = {STARPU_R, STARPU_W},
342 : : .name = "twist1_1d_codelet"
343 : : };
344 : :
345 : : static struct starpu_codelet STARPUFFT(fft1_1d_codelet) = {
346 : : .where =
347 : : #ifdef __STARPU_USE_CUDA
348 : : STARPU_CUDA|
349 : : #endif
350 : : #ifdef STARPU_HAVE_FFTW
351 : : STARPU_CPU|
352 : : #endif
353 : : 0,
354 : : #ifdef __STARPU_USE_CUDA
355 : : .cuda_funcs = {STARPUFFT(fft1_1d_kernel_gpu), NULL},
356 : : #endif
357 : : #ifdef STARPU_HAVE_FFTW
358 : : .cpu_funcs = {STARPUFFT(fft1_1d_kernel_cpu), NULL},
359 : : #endif
360 : : CAN_EXECUTE
361 : : .model = &STARPUFFT(fft1_1d_model),
362 : : .nbuffers = 3,
363 : : .modes = {STARPU_R, STARPU_W, STARPU_R},
364 : : .name = "fft1_1d_codelet"
365 : : };
366 : :
367 : : static struct starpu_codelet STARPUFFT(twist2_1d_codelet) = {
368 : : .where = STARPU_CPU,
369 : : .cpu_funcs = {STARPUFFT(twist2_1d_kernel_cpu), NULL},
370 : : CAN_EXECUTE
371 : : .model = &STARPUFFT(twist2_1d_model),
372 : : .nbuffers = 1,
373 : : .modes = {STARPU_W},
374 : : .name = "twist2_1d_codelet"
375 : : };
376 : :
377 : : static struct starpu_codelet STARPUFFT(fft2_1d_codelet) = {
378 : : .where =
379 : : #ifdef __STARPU_USE_CUDA
380 : : STARPU_CUDA|
381 : : #endif
382 : : #ifdef STARPU_HAVE_FFTW
383 : : STARPU_CPU|
384 : : #endif
385 : : 0,
386 : : #ifdef __STARPU_USE_CUDA
387 : : .cuda_funcs = {STARPUFFT(fft2_1d_kernel_gpu), NULL},
388 : : #endif
389 : : #ifdef STARPU_HAVE_FFTW
390 : : .cpu_funcs = {STARPUFFT(fft2_1d_kernel_cpu), NULL},
391 : : #endif
392 : : CAN_EXECUTE
393 : : .model = &STARPUFFT(fft2_1d_model),
394 : : .nbuffers = 2,
395 : : .modes = {STARPU_R, STARPU_W},
396 : : .name = "fft2_1d_codelet"
397 : : };
398 : :
399 : : static struct starpu_codelet STARPUFFT(twist3_1d_codelet) = {
400 : : .where = STARPU_CPU,
401 : : .cpu_funcs = {STARPUFFT(twist3_1d_kernel_cpu), NULL},
402 : : CAN_EXECUTE
403 : : .model = &STARPUFFT(twist3_1d_model),
404 : : .nbuffers = 1,
405 : : .modes = {STARPU_R},
406 : : .name = "twist3_1d_codelet"
407 : : };
408 : :
409 : : /*
410 : : *
411 : : * Sequential version
412 : : *
413 : : */
414 : :
415 : : #ifdef __STARPU_USE_CUDA
416 : : /* Perform one fft of size n */
417 : : static void
418 : 6 : STARPUFFT(fft_1d_plan_gpu)(void *args)
419 : : {
420 : 6 : STARPUFFT(plan) plan = args;
421 : : cufftResult cures;
422 : 6 : int n = plan->n[0];
423 : 6 : int workerid = starpu_worker_get_id();
424 : :
425 : 6 : cures = cufftPlan1d(&plan->plans[workerid].plan_cuda, n, _CUFFT_C2C, 1);
426 [ - + ]: 6 : if (cures != CUFFT_SUCCESS)
427 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
428 : 6 : cufftSetStream(plan->plans[workerid].plan_cuda, starpu_cuda_get_local_stream());
429 [ - + ]: 6 : if (cures != CUFFT_SUCCESS)
430 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
431 : 6 : }
432 : :
433 : : static void
434 : 2 : STARPUFFT(fft_1d_kernel_gpu)(void *descr[], void *args)
435 : : {
436 : 2 : STARPUFFT(plan) plan = args;
437 : : cufftResult cures;
438 : :
439 : 2 : _cufftComplex * restrict in = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[0]);
440 : 2 : _cufftComplex * restrict out = (_cufftComplex *)STARPU_VECTOR_GET_PTR(descr[1]);
441 : :
442 : 2 : int workerid = starpu_worker_get_id();
443 : :
444 : 2 : task_per_worker[workerid]++;
445 : :
446 [ + - ]: 2 : cures = _cufftExecC2C(plan->plans[workerid].plan_cuda, in, out, plan->sign == -1 ? CUFFT_FORWARD : CUFFT_INVERSE);
447 [ - + ]: 2 : if (cures != CUFFT_SUCCESS)
448 : 0 : STARPU_CUFFT_REPORT_ERROR(cures);
449 : :
450 : 2 : cudaStreamSynchronize(starpu_cuda_get_local_stream());
451 : 2 : }
452 : : #endif
453 : :
454 : : #ifdef STARPU_HAVE_FFTW
455 : : /* Perform one fft of size n */
456 : : static void
457 : 2 : STARPUFFT(fft_1d_kernel_cpu)(void *descr[], void *_args)
458 : : {
459 : 2 : STARPUFFT(plan) plan = _args;
460 : 2 : int workerid = starpu_worker_get_id();
461 : :
462 : 2 : task_per_worker[workerid]++;
463 : :
464 : 2 : STARPUFFT(complex) * restrict in = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[0]);
465 : 2 : STARPUFFT(complex) * restrict out = (STARPUFFT(complex) *)STARPU_VECTOR_GET_PTR(descr[1]);
466 : :
467 : 2 : _FFTW(execute_dft)(plan->plans[workerid].plan_cpu, in, out);
468 : 2 : }
469 : : #endif
470 : :
471 : : static struct starpu_perfmodel STARPUFFT(fft_1d_model) = {
472 : : .type = STARPU_HISTORY_BASED,
473 : : .symbol = TYPE"fft_1d"
474 : : };
475 : :
476 : : static struct starpu_codelet STARPUFFT(fft_1d_codelet) = {
477 : : .where =
478 : : #ifdef __STARPU_USE_CUDA
479 : : STARPU_CUDA|
480 : : #endif
481 : : #ifdef STARPU_HAVE_FFTW
482 : : STARPU_CPU|
483 : : #endif
484 : : 0,
485 : : #ifdef __STARPU_USE_CUDA
486 : : .cuda_funcs = {STARPUFFT(fft_1d_kernel_gpu), NULL},
487 : : #endif
488 : : #ifdef STARPU_HAVE_FFTW
489 : : .cpu_funcs = {STARPUFFT(fft_1d_kernel_cpu), NULL},
490 : : #endif
491 : : CAN_EXECUTE
492 : : .model = &STARPUFFT(fft_1d_model),
493 : : .nbuffers = 2,
494 : : .modes = {STARPU_R, STARPU_W},
495 : : .name = "fft_1d_codelet"
496 : : };
497 : :
498 : : /* Planning:
499 : : *
500 : : * - For each CPU worker, we need to plan the two fftw stages.
501 : : * - For GPU workers, we need to do the planning in the CUDA context, so we do
502 : : * this lazily through the initialised1 and initialised2 flags ; TODO: use
503 : : * starpu_execute_on_each_worker instead (done in the omp branch).
504 : : * - We allocate all the temporary buffers and register them to starpu.
505 : : * - We create all the tasks, but do not submit them yet. It will be possible
506 : : * to reuse them at will to perform several ffts with the same planning.
507 : : */
508 : : STARPUFFT(plan)
509 : 2 : STARPUFFT(plan_dft_1d)(int n, int sign, unsigned flags)
510 : : {
511 : : int workerid;
512 : 2 : int n1 = DIV_1D;
513 : 2 : int n2 = n / n1;
514 : : int n3;
515 : : int z;
516 : : struct starpu_task *task;
517 : :
518 : : if (PARALLEL) {
519 : : #ifdef __STARPU_USE_CUDA
520 : : /* cufft 1D limited to 8M elements */
521 : : while (n2 > 8 << 20) {
522 : : n1 *= 2;
523 : : n2 /= 2;
524 : : }
525 : : #endif
526 : : STARPU_ASSERT(n == n1*n2);
527 : : STARPU_ASSERT(n1 < (1ULL << I_BITS));
528 : :
529 : : /* distribute the n2 second ffts into DIV_1D packages */
530 : : n3 = n2 / DIV_1D;
531 : : STARPU_ASSERT(n2 == n3*DIV_1D);
532 : : }
533 : :
534 : : /* TODO: flags? Automatically set FFTW_MEASURE on calibration? */
535 [ - + ]: 2 : STARPU_ASSERT(flags == 0);
536 : :
537 : 2 : STARPUFFT(plan) plan = malloc(sizeof(*plan));
538 : 2 : memset(plan, 0, sizeof(*plan));
539 : :
540 : : if (PARALLEL) {
541 : : plan->number = STARPU_ATOMIC_ADD(&starpufft_last_plan_number, 1) - 1;
542 : :
543 : : /* The plan number has a limited size */
544 : : STARPU_ASSERT(plan->number < (1ULL << NUMBER_BITS));
545 : : }
546 : :
547 : : /* Just one dimension */
548 : 2 : plan->dim = 1;
549 : 2 : plan->n = malloc(plan->dim * sizeof(*plan->n));
550 : 2 : plan->n[0] = n;
551 : :
552 : : if (PARALLEL) {
553 : : check_dims(plan);
554 : :
555 : : plan->n1 = malloc(plan->dim * sizeof(*plan->n1));
556 : : plan->n1[0] = n1;
557 : : plan->n2 = malloc(plan->dim * sizeof(*plan->n2));
558 : : plan->n2[0] = n2;
559 : : }
560 : :
561 : : /* Note: this is for coherency with the 2D case */
562 : 2 : plan->totsize = n;
563 : :
564 : : if (PARALLEL) {
565 : : plan->totsize1 = n1;
566 : : plan->totsize2 = n2;
567 : : plan->totsize3 = DIV_1D;
568 : : plan->totsize4 = plan->totsize / plan->totsize3;
569 : : }
570 : 2 : plan->type = C2C;
571 : 2 : plan->sign = sign;
572 : :
573 : : if (PARALLEL) {
574 : : /* Compute the w^k just once. */
575 : : compute_roots(plan);
576 : : }
577 : :
578 : : /* Initialize per-worker working set */
579 [ + + ]: 18 : for (workerid = 0; workerid < starpu_worker_get_count(); workerid++) {
580 [ + + - ]: 16 : switch (starpu_worker_get_type(workerid)) {
581 : : case STARPU_CPU_WORKER:
582 : : #ifdef STARPU_HAVE_FFTW
583 : : if (PARALLEL) {
584 : : /* first fft plan: one fft of size n2.
585 : : * FFTW imposes that buffer pointers are known at
586 : : * planning time. */
587 : : plan->plans[workerid].plan1_cpu = _FFTW(plan_dft_1d)(n2, NULL, (void*) 1, sign, _FFTW_FLAGS);
588 : : STARPU_ASSERT(plan->plans[workerid].plan1_cpu);
589 : :
590 : : /* second fft plan: n3 ffts of size n1 */
591 : : plan->plans[workerid].plan2_cpu = _FFTW(plan_many_dft)(plan->dim,
592 : : plan->n1, n3,
593 : : NULL, NULL, 1, plan->totsize1,
594 : : (void*) 1, NULL, 1, plan->totsize1,
595 : : sign, _FFTW_FLAGS);
596 : : STARPU_ASSERT(plan->plans[workerid].plan2_cpu);
597 : : } else {
598 : : /* fft plan: one fft of size n. */
599 : 10 : plan->plans[workerid].plan_cpu = _FFTW(plan_dft_1d)(n, NULL, (void*) 1, sign, _FFTW_FLAGS);
600 [ - + ]: 10 : STARPU_ASSERT(plan->plans[workerid].plan_cpu);
601 : : }
602 : : #else
603 : : /* #warning libstarpufft can not work correctly if libfftw3 is not installed */
604 : : #endif
605 : 10 : break;
606 : : case STARPU_CUDA_WORKER:
607 : 6 : break;
608 : : default:
609 : : /* Do not care, we won't be executing anything there. */
610 : 0 : break;
611 : : }
612 : : }
613 : : #ifdef __STARPU_USE_CUDA
614 : : if (PARALLEL) {
615 : : starpu_execute_on_each_worker(STARPUFFT(fft1_1d_plan_gpu), plan, STARPU_CUDA);
616 : : starpu_execute_on_each_worker(STARPUFFT(fft2_1d_plan_gpu), plan, STARPU_CUDA);
617 : : } else {
618 : 2 : starpu_execute_on_each_worker(STARPUFFT(fft_1d_plan_gpu), plan, STARPU_CUDA);
619 : : }
620 : : #endif
621 : :
622 : : if (PARALLEL) {
623 : : /* Allocate buffers. */
624 : : plan->twisted1 = STARPUFFT(malloc)(plan->totsize * sizeof(*plan->twisted1));
625 : : memset(plan->twisted1, 0, plan->totsize * sizeof(*plan->twisted1));
626 : : plan->fft1 = STARPUFFT(malloc)(plan->totsize * sizeof(*plan->fft1));
627 : : memset(plan->fft1, 0, plan->totsize * sizeof(*plan->fft1));
628 : : plan->twisted2 = STARPUFFT(malloc)(plan->totsize * sizeof(*plan->twisted2));
629 : : memset(plan->twisted2, 0, plan->totsize * sizeof(*plan->twisted2));
630 : : plan->fft2 = STARPUFFT(malloc)(plan->totsize * sizeof(*plan->fft2));
631 : : memset(plan->fft2, 0, plan->totsize * sizeof(*plan->fft2));
632 : :
633 : : /* Allocate handle arrays */
634 : : plan->twisted1_handle = malloc(plan->totsize1 * sizeof(*plan->twisted1_handle));
635 : : plan->fft1_handle = malloc(plan->totsize1 * sizeof(*plan->fft1_handle));
636 : : plan->twisted2_handle = malloc(plan->totsize3 * sizeof(*plan->twisted2_handle));
637 : : plan->fft2_handle = malloc(plan->totsize3 * sizeof(*plan->fft2_handle));
638 : :
639 : : /* Allocate task arrays */
640 : : plan->twist1_tasks = malloc(plan->totsize1 * sizeof(*plan->twist1_tasks));
641 : : plan->fft1_tasks = malloc(plan->totsize1 * sizeof(*plan->fft1_tasks));
642 : : plan->twist2_tasks = malloc(plan->totsize3 * sizeof(*plan->twist2_tasks));
643 : : plan->fft2_tasks = malloc(plan->totsize3 * sizeof(*plan->fft2_tasks));
644 : : plan->twist3_tasks = malloc(plan->totsize3 * sizeof(*plan->twist3_tasks));
645 : :
646 : : /* Allocate codelet argument arrays */
647 : : plan->fft1_args = malloc(plan->totsize1 * sizeof(*plan->fft1_args));
648 : : plan->fft2_args = malloc(plan->totsize3 * sizeof(*plan->fft2_args));
649 : :
650 : : /* Create first-round tasks: DIV_1D tasks of type twist1 and fft1 */
651 : : for (z = 0; z < plan->totsize1; z++) {
652 : : int i = z;
653 : : #define STEP_TAG(step) STEP_TAG_1D(plan, step, i)
654 : :
655 : : /* TODO: get rid of tags */
656 : :
657 : : plan->fft1_args[z].plan = plan;
658 : : plan->fft1_args[z].i = i;
659 : :
660 : : /* Register the twisted1 buffer of size n2. */
661 : : starpu_vector_data_register(&plan->twisted1_handle[z], 0, (uintptr_t) &plan->twisted1[z*plan->totsize2], plan->totsize2, sizeof(*plan->twisted1));
662 : : /* Register the fft1 buffer of size n2. */
663 : : starpu_vector_data_register(&plan->fft1_handle[z], 0, (uintptr_t) &plan->fft1[z*plan->totsize2], plan->totsize2, sizeof(*plan->fft1));
664 : :
665 : : /* We'll need the result of fft1 on the CPU for the second
666 : : * twist anyway, so tell starpu to not keep the fft1 buffer in
667 : : * the GPU. */
668 : : starpu_data_set_wt_mask(plan->fft1_handle[z], 1<<0);
669 : :
670 : : /* Create twist1 task */
671 : : plan->twist1_tasks[z] = task = starpu_task_create();
672 : : task->cl = &STARPUFFT(twist1_1d_codelet);
673 : : /* task->handles[0] = to be filled at execution to point
674 : : to the application input. */
675 : : task->handles[1] = plan->twisted1_handle[z];
676 : : task->cl_arg = &plan->fft1_args[z];
677 : : task->tag_id = STEP_TAG(TWIST1);
678 : : task->use_tag = 1;
679 : : task->destroy = 0;
680 : :
681 : : /* Tell that fft1 depends on twisted1 */
682 : : starpu_tag_declare_deps(STEP_TAG(FFT1),
683 : : 1, STEP_TAG(TWIST1));
684 : :
685 : : /* Create FFT1 task */
686 : : plan->fft1_tasks[z] = task = starpu_task_create();
687 : : task->cl = &STARPUFFT(fft1_1d_codelet);
688 : : task->handles[0] = plan->twisted1_handle[z];
689 : : task->handles[1] = plan->fft1_handle[z];
690 : : task->handles[2] = plan->roots_handle[0];
691 : : task->cl_arg = &plan->fft1_args[z];
692 : : task->tag_id = STEP_TAG(FFT1);
693 : : task->use_tag = 1;
694 : : task->destroy = 0;
695 : :
696 : : /* Tell that the join task will depend on the fft1 task. */
697 : : starpu_tag_declare_deps(STEP_TAG_1D(plan, JOIN, 0),
698 : : 1, STEP_TAG(FFT1));
699 : : #undef STEP_TAG
700 : : }
701 : :
702 : : /* Create the join task, only serving as a dependency point between
703 : : * fft1 and twist2 tasks */
704 : : plan->join_task = task = starpu_task_create();
705 : : task->cl = NULL;
706 : : task->tag_id = STEP_TAG_1D(plan, JOIN, 0);
707 : : task->use_tag = 1;
708 : : task->destroy = 0;
709 : :
710 : : /* Create second-round tasks: DIV_1D batches of n2/DIV_1D twist2, fft2,
711 : : * and twist3 */
712 : : for (z = 0; z < plan->totsize3; z++) {
713 : : int jj = z;
714 : : #define STEP_TAG(step) STEP_TAG_1D(plan, step, jj)
715 : :
716 : : plan->fft2_args[z].plan = plan;
717 : : plan->fft2_args[z].jj = jj;
718 : :
719 : : /* Register n3 twisted2 buffers of size n1 */
720 : : starpu_vector_data_register(&plan->twisted2_handle[z], 0, (uintptr_t) &plan->twisted2[z*plan->totsize4], plan->totsize4, sizeof(*plan->twisted2));
721 : : starpu_vector_data_register(&plan->fft2_handle[z], 0, (uintptr_t) &plan->fft2[z*plan->totsize4], plan->totsize4, sizeof(*plan->fft2));
722 : :
723 : : /* We'll need the result of fft2 on the CPU for the third
724 : : * twist anyway, so tell starpu to not keep the fft2 buffer in
725 : : * the GPU. */
726 : : starpu_data_set_wt_mask(plan->fft2_handle[z], 1<<0);
727 : :
728 : : /* Tell that twisted2 depends on the join task */
729 : : starpu_tag_declare_deps(STEP_TAG(TWIST2),
730 : : 1, STEP_TAG_1D(plan, JOIN, 0));
731 : :
732 : : /* Create twist2 task */
733 : : plan->twist2_tasks[z] = task = starpu_task_create();
734 : : task->cl = &STARPUFFT(twist2_1d_codelet);
735 : : task->handles[0] = plan->twisted2_handle[z];
736 : : task->cl_arg = &plan->fft2_args[z];
737 : : task->tag_id = STEP_TAG(TWIST2);
738 : : task->use_tag = 1;
739 : : task->destroy = 0;
740 : :
741 : : /* Tell that fft2 depends on twisted2 */
742 : : starpu_tag_declare_deps(STEP_TAG(FFT2),
743 : : 1, STEP_TAG(TWIST2));
744 : :
745 : : /* Create FFT2 task */
746 : : plan->fft2_tasks[z] = task = starpu_task_create();
747 : : task->cl = &STARPUFFT(fft2_1d_codelet);
748 : : task->handles[0] = plan->twisted2_handle[z];
749 : : task->handles[1] = plan->fft2_handle[z];
750 : : task->cl_arg = &plan->fft2_args[z];
751 : : task->tag_id = STEP_TAG(FFT2);
752 : : task->use_tag = 1;
753 : : task->destroy = 0;
754 : :
755 : : /* Tell that twist3 depends on fft2 */
756 : : starpu_tag_declare_deps(STEP_TAG(TWIST3),
757 : : 1, STEP_TAG(FFT2));
758 : :
759 : : /* Create twist3 tasks */
760 : : /* These run only on CPUs and thus write directly into the
761 : : * application output buffer. */
762 : : plan->twist3_tasks[z] = task = starpu_task_create();
763 : : task->cl = &STARPUFFT(twist3_1d_codelet);
764 : : task->handles[0] = plan->fft2_handle[z];
765 : : task->cl_arg = &plan->fft2_args[z];
766 : : task->tag_id = STEP_TAG(TWIST3);
767 : : task->use_tag = 1;
768 : : task->destroy = 0;
769 : :
770 : : /* Tell that to be completely finished we need to have finished
771 : : * this twisted3 */
772 : : starpu_tag_declare_deps(STEP_TAG_1D(plan, END, 0),
773 : : 1, STEP_TAG(TWIST3));
774 : : #undef STEP_TAG
775 : : }
776 : :
777 : : /* Create end task, only serving as a join point. */
778 : : plan->end_task = task = starpu_task_create();
779 : : task->cl = NULL;
780 : : task->tag_id = STEP_TAG_1D(plan, END, 0);
781 : : task->use_tag = 1;
782 : : task->destroy = 0;
783 : :
784 : : }
785 : :
786 : 2 : return plan;
787 : : }
788 : :
789 : : /* Actually submit all the tasks. */
790 : : static struct starpu_task *
791 : 4 : STARPUFFT(start1dC2C)(STARPUFFT(plan) plan, starpu_data_handle_t in, starpu_data_handle_t out)
792 : : {
793 [ - + ]: 4 : STARPU_ASSERT(plan->type == C2C);
794 : : int z;
795 : : int ret;
796 : :
797 : : if (PARALLEL) {
798 : : for (z=0; z < plan->totsize1; z++) {
799 : : ret = starpu_task_submit(plan->twist1_tasks[z]);
800 : : if (ret == -ENODEV) return NULL;
801 : : STARPU_CHECK_RETURN_VALUE(ret, "starpu_task_submit");
802 : : ret = starpu_task_submit(plan->fft1_tasks[z]);
803 : : if (ret == -ENODEV) return NULL;
804 : : STARPU_CHECK_RETURN_VALUE(ret, "starpu_task_submit");
805 : : }
806 : :
807 : : ret = starpu_task_submit(plan->join_task);
808 : : if (ret == -ENODEV) return NULL;
809 : : STARPU_CHECK_RETURN_VALUE(ret, "starpu_task_submit");
810 : :
811 : : for (z=0; z < plan->totsize3; z++) {
812 : : ret = starpu_task_submit(plan->twist2_tasks[z]);
813 : : if (ret == -ENODEV) return NULL;
814 : : STARPU_CHECK_RETURN_VALUE(ret, "starpu_task_submit");
815 : : ret = starpu_task_submit(plan->fft2_tasks[z]);
816 : : if (ret == -ENODEV) return NULL;
817 : : STARPU_CHECK_RETURN_VALUE(ret, "starpu_task_submit");
818 : : ret = starpu_task_submit(plan->twist3_tasks[z]);
819 : : if (ret == -ENODEV) return NULL;
820 : : STARPU_CHECK_RETURN_VALUE(ret, "starpu_task_submit");
821 : : }
822 : :
823 : : ret = starpu_task_submit(plan->end_task);
824 : : if (ret == -ENODEV) return NULL;
825 : : STARPU_CHECK_RETURN_VALUE(ret, "starpu_task_submit");
826 : :
827 : : return plan->end_task;
828 : : } else /* !PARALLEL */ {
829 : : struct starpu_task *task;
830 : :
831 : : /* Create FFT task */
832 : 4 : task = starpu_task_create();
833 : 4 : task->detach = 0;
834 : 4 : task->cl = &STARPUFFT(fft_1d_codelet);
835 : 4 : task->handles[0] = in;
836 : 4 : task->handles[1] = out;
837 : 4 : task->cl_arg = plan;
838 : :
839 : 4 : ret = starpu_task_submit(task);
840 [ - + ]: 4 : if (ret == -ENODEV) return NULL;
841 [ - + ]: 4 : STARPU_CHECK_RETURN_VALUE(ret, "starpu_task_submit");
842 : 4 : return task;
843 : : }
844 : : }
845 : :
846 : : /* Free all the tags. The generic code handles freeing the buffers. */
847 : : static void
848 : 0 : STARPUFFT(free_1d_tags)(STARPUFFT(plan) plan)
849 : : {
850 : : unsigned i;
851 : 0 : int n1 = plan->n1[0];
852 : :
853 : : if (!PARALLEL)
854 : 0 : return;
855 : :
856 : : for (i = 0; i < n1; i++) {
857 : : starpu_tag_remove(STEP_TAG_1D(plan, TWIST1, i));
858 : : starpu_tag_remove(STEP_TAG_1D(plan, FFT1, i));
859 : : }
860 : :
861 : : starpu_tag_remove(STEP_TAG_1D(plan, JOIN, 0));
862 : :
863 : : for (i = 0; i < DIV_1D; i++) {
864 : : starpu_tag_remove(STEP_TAG_1D(plan, TWIST2, i));
865 : : starpu_tag_remove(STEP_TAG_1D(plan, FFT2, i));
866 : : starpu_tag_remove(STEP_TAG_1D(plan, TWIST3, i));
867 : : }
868 : :
869 : : starpu_tag_remove(STEP_TAG_1D(plan, END, 0));
870 : : }
|
