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
22 #include "dft/ct.h"
23
24 ct_solver *(*X(mksolver_ct_hook))(size_t, INT, int,
25 ct_mkinferior, ct_force_vrecursion) = 0;
26
27 typedef struct {
28 plan_dft super;
29 plan *cld;
30 plan *cldw;
31 INT r;
32 } P;
33
apply_dit(const plan * ego_,R * ri,R * ii,R * ro,R * io)34 static void apply_dit(const plan *ego_, R *ri, R *ii, R *ro, R *io)
35 {
36 const P *ego = (const P *) ego_;
37 plan_dft *cld;
38 plan_dftw *cldw;
39
40 cld = (plan_dft *) ego->cld;
41 cld->apply(ego->cld, ri, ii, ro, io);
42
43 cldw = (plan_dftw *) ego->cldw;
44 cldw->apply(ego->cldw, ro, io);
45 }
46
apply_dif(const plan * ego_,R * ri,R * ii,R * ro,R * io)47 static void apply_dif(const plan *ego_, R *ri, R *ii, R *ro, R *io)
48 {
49 const P *ego = (const P *) ego_;
50 plan_dft *cld;
51 plan_dftw *cldw;
52
53 cldw = (plan_dftw *) ego->cldw;
54 cldw->apply(ego->cldw, ri, ii);
55
56 cld = (plan_dft *) ego->cld;
57 cld->apply(ego->cld, ri, ii, ro, io);
58 }
59
awake(plan * ego_,enum wakefulness wakefulness)60 static void awake(plan *ego_, enum wakefulness wakefulness)
61 {
62 P *ego = (P *) ego_;
63 X(plan_awake)(ego->cld, wakefulness);
64 X(plan_awake)(ego->cldw, wakefulness);
65 }
66
destroy(plan * ego_)67 static void destroy(plan *ego_)
68 {
69 P *ego = (P *) ego_;
70 X(plan_destroy_internal)(ego->cldw);
71 X(plan_destroy_internal)(ego->cld);
72 }
73
print(const plan * ego_,printer * p)74 static void print(const plan *ego_, printer *p)
75 {
76 const P *ego = (const P *) ego_;
77 p->print(p, "(dft-ct-%s/%D%(%p%)%(%p%))",
78 ego->super.apply == apply_dit ? "dit" : "dif",
79 ego->r, ego->cldw, ego->cld);
80 }
81
applicable0(const ct_solver * ego,const problem * p_,planner * plnr)82 static int applicable0(const ct_solver *ego, const problem *p_, planner *plnr)
83 {
84 const problem_dft *p = (const problem_dft *) p_;
85 INT r;
86
87 return (1
88 && p->sz->rnk == 1
89 && p->vecsz->rnk <= 1
90
91 /* DIF destroys the input and we don't like it */
92 && (ego->dec == DECDIT ||
93 p->ri == p->ro ||
94 !NO_DESTROY_INPUTP(plnr))
95
96 && ((r = X(choose_radix)(ego->r, p->sz->dims[0].n)) > 1)
97 && p->sz->dims[0].n > r);
98 }
99
100
X(ct_applicable)101 int X(ct_applicable)(const ct_solver *ego, const problem *p_, planner *plnr)
102 {
103 const problem_dft *p;
104
105 if (!applicable0(ego, p_, plnr))
106 return 0;
107
108 p = (const problem_dft *) p_;
109
110 return (0
111 || ego->dec == DECDIF+TRANSPOSE
112 || p->vecsz->rnk == 0
113 || !NO_VRECURSEP(plnr)
114 || (ego->force_vrecursionp && ego->force_vrecursionp(ego, p))
115 );
116 }
117
118
mkplan(const solver * ego_,const problem * p_,planner * plnr)119 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
120 {
121 const ct_solver *ego = (const ct_solver *) ego_;
122 const problem_dft *p;
123 P *pln = 0;
124 plan *cld = 0, *cldw = 0;
125 INT n, r, m, v, ivs, ovs;
126 iodim *d;
127
128 static const plan_adt padt = {
129 X(dft_solve), awake, print, destroy
130 };
131
132 if ((NO_NONTHREADEDP(plnr)) || !X(ct_applicable)(ego, p_, plnr))
133 return (plan *) 0;
134
135 p = (const problem_dft *) p_;
136 d = p->sz->dims;
137 n = d[0].n;
138 r = X(choose_radix)(ego->r, n);
139 m = n / r;
140
141 X(tensor_tornk1)(p->vecsz, &v, &ivs, &ovs);
142
143 switch (ego->dec) {
144 case DECDIT:
145 {
146 cldw = ego->mkcldw(ego,
147 r, m * d[0].os, m * d[0].os,
148 m, d[0].os,
149 v, ovs, ovs,
150 0, m,
151 p->ro, p->io, plnr);
152 if (!cldw) goto nada;
153
154 cld = X(mkplan_d)(plnr,
155 X(mkproblem_dft_d)(
156 X(mktensor_1d)(m, r * d[0].is, d[0].os),
157 X(mktensor_2d)(r, d[0].is, m * d[0].os,
158 v, ivs, ovs),
159 p->ri, p->ii, p->ro, p->io)
160 );
161 if (!cld) goto nada;
162
163 pln = MKPLAN_DFT(P, &padt, apply_dit);
164 break;
165 }
166 case DECDIF:
167 case DECDIF+TRANSPOSE:
168 {
169 INT cors, covs; /* cldw ors, ovs */
170 if (ego->dec == DECDIF+TRANSPOSE) {
171 cors = ivs;
172 covs = m * d[0].is;
173 /* ensure that we generate well-formed dftw subproblems */
174 /* FIXME: too conservative */
175 if (!(1
176 && r == v
177 && d[0].is == r * cors))
178 goto nada;
179
180 /* FIXME: allow in-place only for now, like in
181 fftw-3.[01] */
182 if (!(1
183 && p->ri == p->ro
184 && d[0].is == r * d[0].os
185 && cors == d[0].os
186 && covs == ovs
187 ))
188 goto nada;
189 } else {
190 cors = m * d[0].is;
191 covs = ivs;
192 }
193
194 cldw = ego->mkcldw(ego,
195 r, m * d[0].is, cors,
196 m, d[0].is,
197 v, ivs, covs,
198 0, m,
199 p->ri, p->ii, plnr);
200 if (!cldw) goto nada;
201
202 cld = X(mkplan_d)(plnr,
203 X(mkproblem_dft_d)(
204 X(mktensor_1d)(m, d[0].is, r * d[0].os),
205 X(mktensor_2d)(r, cors, d[0].os,
206 v, covs, ovs),
207 p->ri, p->ii, p->ro, p->io)
208 );
209 if (!cld) goto nada;
210
211 pln = MKPLAN_DFT(P, &padt, apply_dif);
212 break;
213 }
214
215 default: A(0);
216
217 }
218
219 pln->cld = cld;
220 pln->cldw = cldw;
221 pln->r = r;
222 X(ops_add)(&cld->ops, &cldw->ops, &pln->super.super.ops);
223
224 /* inherit could_prune_now_p attribute from cldw */
225 pln->super.super.could_prune_now_p = cldw->could_prune_now_p;
226 return &(pln->super.super);
227
228 nada:
229 X(plan_destroy_internal)(cldw);
230 X(plan_destroy_internal)(cld);
231 return (plan *) 0;
232 }
233
X(mksolver_ct)234 ct_solver *X(mksolver_ct)(size_t size, INT r, int dec,
235 ct_mkinferior mkcldw,
236 ct_force_vrecursion force_vrecursionp)
237 {
238 static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 };
239 ct_solver *slv = (ct_solver *)X(mksolver)(size, &sadt);
240 slv->r = r;
241 slv->dec = dec;
242 slv->mkcldw = mkcldw;
243 slv->force_vrecursionp = force_vrecursionp;
244 return slv;
245 }
246
X(mkplan_dftw)247 plan *X(mkplan_dftw)(size_t size, const plan_adt *adt, dftwapply apply)
248 {
249 plan_dftw *ego;
250
251 ego = (plan_dftw *) X(mkplan)(size, adt);
252 ego->apply = apply;
253
254 return &(ego->super);
255 }
256