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 "rdft/hc2hc.h"
22
23 hc2hc_solver *(*X(mksolver_hc2hc_hook))(size_t, INT, hc2hc_mkinferior) = 0;
24
25 typedef struct {
26 plan_rdft super;
27 plan *cld;
28 plan *cldw;
29 INT r;
30 } P;
31
apply_dit(const plan * ego_,R * I,R * O)32 static void apply_dit(const plan *ego_, R *I, R *O)
33 {
34 const P *ego = (const P *) ego_;
35 plan_rdft *cld;
36 plan_hc2hc *cldw;
37
38 cld = (plan_rdft *) ego->cld;
39 cld->apply(ego->cld, I, O);
40
41 cldw = (plan_hc2hc *) ego->cldw;
42 cldw->apply(ego->cldw, O);
43 }
44
apply_dif(const plan * ego_,R * I,R * O)45 static void apply_dif(const plan *ego_, R *I, R *O)
46 {
47 const P *ego = (const P *) ego_;
48 plan_rdft *cld;
49 plan_hc2hc *cldw;
50
51 cldw = (plan_hc2hc *) ego->cldw;
52 cldw->apply(ego->cldw, I);
53
54 cld = (plan_rdft *) ego->cld;
55 cld->apply(ego->cld, I, O);
56 }
57
awake(plan * ego_,enum wakefulness wakefulness)58 static void awake(plan *ego_, enum wakefulness wakefulness)
59 {
60 P *ego = (P *) ego_;
61 X(plan_awake)(ego->cld, wakefulness);
62 X(plan_awake)(ego->cldw, wakefulness);
63 }
64
destroy(plan * ego_)65 static void destroy(plan *ego_)
66 {
67 P *ego = (P *) ego_;
68 X(plan_destroy_internal)(ego->cldw);
69 X(plan_destroy_internal)(ego->cld);
70 }
71
print(const plan * ego_,printer * p)72 static void print(const plan *ego_, printer *p)
73 {
74 const P *ego = (const P *) ego_;
75 p->print(p, "(rdft-ct-%s/%D%(%p%)%(%p%))",
76 ego->super.apply == apply_dit ? "dit" : "dif",
77 ego->r, ego->cldw, ego->cld);
78 }
79
applicable0(const hc2hc_solver * ego,const problem * p_,planner * plnr)80 static int applicable0(const hc2hc_solver *ego, const problem *p_, planner *plnr)
81 {
82 const problem_rdft *p = (const problem_rdft *) p_;
83 INT r;
84
85 return (1
86 && p->sz->rnk == 1
87 && p->vecsz->rnk <= 1
88
89 && (/* either the problem is R2HC, which is solved by DIT */
90 (p->kind[0] == R2HC)
91 ||
92 /* or the problem is HC2R, in which case it is solved
93 by DIF, which destroys the input */
94 (p->kind[0] == HC2R &&
95 (p->I == p->O || !NO_DESTROY_INPUTP(plnr))))
96
97 && ((r = X(choose_radix)(ego->r, p->sz->dims[0].n)) > 0)
98 && p->sz->dims[0].n > r);
99 }
100
X(hc2hc_applicable)101 int X(hc2hc_applicable)(const hc2hc_solver *ego, const problem *p_, planner *plnr)
102 {
103 const problem_rdft *p;
104
105 if (!applicable0(ego, p_, plnr))
106 return 0;
107
108 p = (const problem_rdft *) p_;
109
110 return (0
111 || p->vecsz->rnk == 0
112 || !NO_VRECURSEP(plnr)
113 );
114 }
115
mkplan(const solver * ego_,const problem * p_,planner * plnr)116 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
117 {
118 const hc2hc_solver *ego = (const hc2hc_solver *) ego_;
119 const problem_rdft *p;
120 P *pln = 0;
121 plan *cld = 0, *cldw = 0;
122 INT n, r, m, v, ivs, ovs;
123 iodim *d;
124
125 static const plan_adt padt = {
126 X(rdft_solve), awake, print, destroy
127 };
128
129 if (NO_NONTHREADEDP(plnr) || !X(hc2hc_applicable)(ego, p_, plnr))
130 return (plan *) 0;
131
132 p = (const problem_rdft *) p_;
133 d = p->sz->dims;
134 n = d[0].n;
135 r = X(choose_radix)(ego->r, n);
136 m = n / r;
137
138 X(tensor_tornk1)(p->vecsz, &v, &ivs, &ovs);
139
140 switch (p->kind[0]) {
141 case R2HC:
142 cldw = ego->mkcldw(ego,
143 R2HC, r, m, d[0].os, v, ovs, 0, (m+2)/2,
144 p->O, plnr);
145 if (!cldw) goto nada;
146
147 cld = X(mkplan_d)(plnr,
148 X(mkproblem_rdft_d)(
149 X(mktensor_1d)(m, r * d[0].is, d[0].os),
150 X(mktensor_2d)(r, d[0].is, m * d[0].os,
151 v, ivs, ovs),
152 p->I, p->O, p->kind)
153 );
154 if (!cld) goto nada;
155
156 pln = MKPLAN_RDFT(P, &padt, apply_dit);
157 break;
158
159 case HC2R:
160 cldw = ego->mkcldw(ego,
161 HC2R, r, m, d[0].is, v, ivs, 0, (m+2)/2,
162 p->I, plnr);
163 if (!cldw) goto nada;
164
165 cld = X(mkplan_d)(plnr,
166 X(mkproblem_rdft_d)(
167 X(mktensor_1d)(m, d[0].is, r * d[0].os),
168 X(mktensor_2d)(r, m * d[0].is, d[0].os,
169 v, ivs, ovs),
170 p->I, p->O, p->kind)
171 );
172 if (!cld) goto nada;
173
174 pln = MKPLAN_RDFT(P, &padt, apply_dif);
175 break;
176
177 default:
178 A(0);
179 }
180
181 pln->cld = cld;
182 pln->cldw = cldw;
183 pln->r = r;
184 X(ops_add)(&cld->ops, &cldw->ops, &pln->super.super.ops);
185
186 /* inherit could_prune_now_p attribute from cldw */
187 pln->super.super.could_prune_now_p = cldw->could_prune_now_p;
188
189 return &(pln->super.super);
190
191 nada:
192 X(plan_destroy_internal)(cldw);
193 X(plan_destroy_internal)(cld);
194 return (plan *) 0;
195 }
196
X(mksolver_hc2hc)197 hc2hc_solver *X(mksolver_hc2hc)(size_t size, INT r, hc2hc_mkinferior mkcldw)
198 {
199 static const solver_adt sadt = { PROBLEM_RDFT, mkplan, 0 };
200 hc2hc_solver *slv = (hc2hc_solver *)X(mksolver)(size, &sadt);
201 slv->r = r;
202 slv->mkcldw = mkcldw;
203 return slv;
204 }
205
X(mkplan_hc2hc)206 plan *X(mkplan_hc2hc)(size_t size, const plan_adt *adt, hc2hcapply apply)
207 {
208 plan_hc2hc *ego;
209
210 ego = (plan_hc2hc *) X(mkplan)(size, adt);
211 ego->apply = apply;
212
213 return &(ego->super);
214 }
215