1 /*
2 * Copyright (c) 2003, 2007-14 Matteo Frigo
3 * Copyright (c) 2003, 2007-14 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 *
19 */
20
21 #include "dft/dft.h"
22
23 typedef struct {
24 solver super;
25 } S;
26
27 typedef struct {
28 plan_dft super;
29 INT n; /* problem size */
30 INT nb; /* size of convolution */
31 R *w; /* lambda k . exp(2*pi*i*k^2/(2*n)) */
32 R *W; /* DFT(w) */
33 plan *cldf;
34 INT is, os;
35 } P;
36
bluestein_sequence(enum wakefulness wakefulness,INT n,R * w)37 static void bluestein_sequence(enum wakefulness wakefulness, INT n, R *w)
38 {
39 INT k, ksq, n2 = 2 * n;
40 triggen *t = X(mktriggen)(wakefulness, n2);
41
42 ksq = 0;
43 for (k = 0; k < n; ++k) {
44 t->cexp(t, ksq, w+2*k);
45 /* careful with overflow */
46 ksq += 2*k + 1; while (ksq > n2) ksq -= n2;
47 }
48
49 X(triggen_destroy)(t);
50 }
51
mktwiddle(enum wakefulness wakefulness,P * p)52 static void mktwiddle(enum wakefulness wakefulness, P *p)
53 {
54 INT i;
55 INT n = p->n, nb = p->nb;
56 R *w, *W;
57 E nbf = (E)nb;
58
59 p->w = w = (R *) MALLOC(2 * n * sizeof(R), TWIDDLES);
60 p->W = W = (R *) MALLOC(2 * nb * sizeof(R), TWIDDLES);
61
62 bluestein_sequence(wakefulness, n, w);
63
64 for (i = 0; i < nb; ++i)
65 W[2*i] = W[2*i+1] = K(0.0);
66
67 W[0] = w[0] / nbf;
68 W[1] = w[1] / nbf;
69
70 for (i = 1; i < n; ++i) {
71 W[2*i] = W[2*(nb-i)] = w[2*i] / nbf;
72 W[2*i+1] = W[2*(nb-i)+1] = w[2*i+1] / nbf;
73 }
74
75 {
76 plan_dft *cldf = (plan_dft *)p->cldf;
77 /* cldf must be awake */
78 cldf->apply(p->cldf, W, W+1, W, W+1);
79 }
80 }
81
apply(const plan * ego_,R * ri,R * ii,R * ro,R * io)82 static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io)
83 {
84 const P *ego = (const P *) ego_;
85 INT i, n = ego->n, nb = ego->nb, is = ego->is, os = ego->os;
86 R *w = ego->w, *W = ego->W;
87 R *b = (R *) MALLOC(2 * nb * sizeof(R), BUFFERS);
88
89 /* multiply input by conjugate bluestein sequence */
90 for (i = 0; i < n; ++i) {
91 E xr = ri[i*is], xi = ii[i*is];
92 E wr = w[2*i], wi = w[2*i+1];
93 b[2*i] = xr * wr + xi * wi;
94 b[2*i+1] = xi * wr - xr * wi;
95 }
96
97 for (; i < nb; ++i) b[2*i] = b[2*i+1] = K(0.0);
98
99 /* convolution: FFT */
100 {
101 plan_dft *cldf = (plan_dft *)ego->cldf;
102 cldf->apply(ego->cldf, b, b+1, b, b+1);
103 }
104
105 /* convolution: pointwise multiplication */
106 for (i = 0; i < nb; ++i) {
107 E xr = b[2*i], xi = b[2*i+1];
108 E wr = W[2*i], wi = W[2*i+1];
109 b[2*i] = xi * wr + xr * wi;
110 b[2*i+1] = xr * wr - xi * wi;
111 }
112
113 /* convolution: IFFT by FFT with real/imag input/output swapped */
114 {
115 plan_dft *cldf = (plan_dft *)ego->cldf;
116 cldf->apply(ego->cldf, b, b+1, b, b+1);
117 }
118
119 /* multiply output by conjugate bluestein sequence */
120 for (i = 0; i < n; ++i) {
121 E xi = b[2*i], xr = b[2*i+1];
122 E wr = w[2*i], wi = w[2*i+1];
123 ro[i*os] = xr * wr + xi * wi;
124 io[i*os] = xi * wr - xr * wi;
125 }
126
127 X(ifree)(b);
128 }
129
awake(plan * ego_,enum wakefulness wakefulness)130 static void awake(plan *ego_, enum wakefulness wakefulness)
131 {
132 P *ego = (P *) ego_;
133
134 X(plan_awake)(ego->cldf, wakefulness);
135
136 switch (wakefulness) {
137 case SLEEPY:
138 X(ifree0)(ego->w); ego->w = 0;
139 X(ifree0)(ego->W); ego->W = 0;
140 break;
141 default:
142 A(!ego->w);
143 mktwiddle(wakefulness, ego);
144 break;
145 }
146 }
147
applicable(const solver * ego,const problem * p_,const planner * plnr)148 static int applicable(const solver *ego, const problem *p_,
149 const planner *plnr)
150 {
151 const problem_dft *p = (const problem_dft *) p_;
152 UNUSED(ego);
153 return (1
154 && p->sz->rnk == 1
155 && p->vecsz->rnk == 0
156 /* FIXME: allow other sizes */
157 && X(is_prime)(p->sz->dims[0].n)
158
159 /* FIXME: avoid infinite recursion of bluestein with itself.
160 This works because all factors in child problems are 2, 3, 5 */
161 && p->sz->dims[0].n > 16
162
163 && CIMPLIES(NO_SLOWP(plnr), p->sz->dims[0].n > BLUESTEIN_MAX_SLOW)
164 );
165 }
166
destroy(plan * ego_)167 static void destroy(plan *ego_)
168 {
169 P *ego = (P *) ego_;
170 X(plan_destroy_internal)(ego->cldf);
171 }
172
print(const plan * ego_,printer * p)173 static void print(const plan *ego_, printer *p)
174 {
175 const P *ego = (const P *)ego_;
176 p->print(p, "(dft-bluestein-%D/%D%(%p%))",
177 ego->n, ego->nb, ego->cldf);
178 }
179
choose_transform_size(INT minsz)180 static INT choose_transform_size(INT minsz)
181 {
182 while (!X(factors_into_small_primes)(minsz))
183 ++minsz;
184 return minsz;
185 }
186
mkplan(const solver * ego,const problem * p_,planner * plnr)187 static plan *mkplan(const solver *ego, const problem *p_, planner *plnr)
188 {
189 const problem_dft *p = (const problem_dft *) p_;
190 P *pln;
191 INT n, nb;
192 plan *cldf = 0;
193 R *buf = (R *) 0;
194
195 static const plan_adt padt = {
196 X(dft_solve), awake, print, destroy
197 };
198
199 if (!applicable(ego, p_, plnr))
200 return (plan *) 0;
201
202 n = p->sz->dims[0].n;
203 nb = choose_transform_size(2 * n - 1);
204 buf = (R *) MALLOC(2 * nb * sizeof(R), BUFFERS);
205
206 cldf = X(mkplan_f_d)(plnr,
207 X(mkproblem_dft_d)(X(mktensor_1d)(nb, 2, 2),
208 X(mktensor_1d)(1, 0, 0),
209 buf, buf+1,
210 buf, buf+1),
211 NO_SLOW, 0, 0);
212 if (!cldf) goto nada;
213
214 X(ifree)(buf);
215
216 pln = MKPLAN_DFT(P, &padt, apply);
217
218 pln->n = n;
219 pln->nb = nb;
220 pln->w = 0;
221 pln->W = 0;
222 pln->cldf = cldf;
223 pln->is = p->sz->dims[0].is;
224 pln->os = p->sz->dims[0].os;
225
226 X(ops_add)(&cldf->ops, &cldf->ops, &pln->super.super.ops);
227 pln->super.super.ops.add += 4 * n + 2 * nb;
228 pln->super.super.ops.mul += 8 * n + 4 * nb;
229 pln->super.super.ops.other += 6 * (n + nb);
230
231 return &(pln->super.super);
232
233 nada:
234 X(ifree0)(buf);
235 X(plan_destroy_internal)(cldf);
236 return (plan *)0;
237 }
238
239
mksolver(void)240 static solver *mksolver(void)
241 {
242 static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 };
243 S *slv = MKSOLVER(S, &sadt);
244 return &(slv->super);
245 }
246
X(dft_bluestein_register)247 void X(dft_bluestein_register)(planner *p)
248 {
249 REGISTER_SOLVER(p, mksolver());
250 }
251