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 /* Compute the complex DFT by combining R2HC RDFTs on the real
23 and imaginary parts. This could be useful for people just wanting
24 to link to the real codelets and not the complex ones. It could
25 also even be faster than the complex algorithms for split (as opposed
26 to interleaved) real/imag complex data. */
27
28 #include "rdft/rdft.h"
29 #include "dft/dft.h"
30
31 typedef struct {
32 solver super;
33 } S;
34
35 typedef struct {
36 plan_dft super;
37 plan *cld;
38 INT ishift, oshift;
39 INT os;
40 INT n;
41 } P;
42
apply(const plan * ego_,R * ri,R * ii,R * ro,R * io)43 static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io)
44 {
45 const P *ego = (const P *) ego_;
46 INT n;
47
48 UNUSED(ii);
49
50 { /* transform vector of real & imag parts: */
51 plan_rdft *cld = (plan_rdft *) ego->cld;
52 cld->apply((plan *) cld, ri + ego->ishift, ro + ego->oshift);
53 }
54
55 n = ego->n;
56 if (n > 1) {
57 INT i, os = ego->os;
58 for (i = 1; i < (n + 1)/2; ++i) {
59 E rop, iop, iom, rom;
60 rop = ro[os * i];
61 iop = io[os * i];
62 rom = ro[os * (n - i)];
63 iom = io[os * (n - i)];
64 ro[os * i] = rop - iom;
65 io[os * i] = iop + rom;
66 ro[os * (n - i)] = rop + iom;
67 io[os * (n - i)] = iop - rom;
68 }
69 }
70 }
71
awake(plan * ego_,enum wakefulness wakefulness)72 static void awake(plan *ego_, enum wakefulness wakefulness)
73 {
74 P *ego = (P *) ego_;
75 X(plan_awake)(ego->cld, wakefulness);
76 }
77
destroy(plan * ego_)78 static void destroy(plan *ego_)
79 {
80 P *ego = (P *) ego_;
81 X(plan_destroy_internal)(ego->cld);
82 }
83
print(const plan * ego_,printer * p)84 static void print(const plan *ego_, printer *p)
85 {
86 const P *ego = (const P *) ego_;
87 p->print(p, "(dft-r2hc-%D%(%p%))", ego->n, ego->cld);
88 }
89
90
applicable0(const problem * p_)91 static int applicable0(const problem *p_)
92 {
93 const problem_dft *p = (const problem_dft *) p_;
94 return ((p->sz->rnk == 1 && p->vecsz->rnk == 0)
95 || (p->sz->rnk == 0 && FINITE_RNK(p->vecsz->rnk))
96 );
97 }
98
splitp(R * r,R * i,INT n,INT s)99 static int splitp(R *r, R *i, INT n, INT s)
100 {
101 return ((r > i ? (r - i) : (i - r)) >= n * (s > 0 ? s : 0-s));
102 }
103
applicable(const problem * p_,const planner * plnr)104 static int applicable(const problem *p_, const planner *plnr)
105 {
106 if (!applicable0(p_)) return 0;
107
108 {
109 const problem_dft *p = (const problem_dft *) p_;
110
111 /* rank-0 problems are always OK */
112 if (p->sz->rnk == 0) return 1;
113
114 /* this solver is ok for split arrays */
115 if (p->sz->rnk == 1 &&
116 splitp(p->ri, p->ii, p->sz->dims[0].n, p->sz->dims[0].is) &&
117 splitp(p->ro, p->io, p->sz->dims[0].n, p->sz->dims[0].os))
118 return 1;
119
120 return !(NO_DFT_R2HCP(plnr));
121 }
122 }
123
mkplan(const solver * ego_,const problem * p_,planner * plnr)124 static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
125 {
126 P *pln;
127 const problem_dft *p;
128 plan *cld;
129 INT ishift = 0, oshift = 0;
130
131 static const plan_adt padt = {
132 X(dft_solve), awake, print, destroy
133 };
134
135 UNUSED(ego_);
136 if (!applicable(p_, plnr))
137 return (plan *)0;
138
139 p = (const problem_dft *) p_;
140
141 {
142 tensor *ri_vec = X(mktensor_1d)(2, p->ii - p->ri, p->io - p->ro);
143 tensor *cld_vec = X(tensor_append)(ri_vec, p->vecsz);
144 int i;
145 for (i = 0; i < cld_vec->rnk; ++i) { /* make all istrides > 0 */
146 if (cld_vec->dims[i].is < 0) {
147 INT nm1 = cld_vec->dims[i].n - 1;
148 ishift -= nm1 * (cld_vec->dims[i].is *= -1);
149 oshift -= nm1 * (cld_vec->dims[i].os *= -1);
150 }
151 }
152 cld = X(mkplan_d)(plnr,
153 X(mkproblem_rdft_1)(p->sz, cld_vec,
154 p->ri + ishift,
155 p->ro + oshift, R2HC));
156 X(tensor_destroy2)(ri_vec, cld_vec);
157 }
158 if (!cld) return (plan *)0;
159
160 pln = MKPLAN_DFT(P, &padt, apply);
161
162 if (p->sz->rnk == 0) {
163 pln->n = 1;
164 pln->os = 0;
165 }
166 else {
167 pln->n = p->sz->dims[0].n;
168 pln->os = p->sz->dims[0].os;
169 }
170 pln->ishift = ishift;
171 pln->oshift = oshift;
172
173 pln->cld = cld;
174
175 pln->super.super.ops = cld->ops;
176 pln->super.super.ops.other += 8 * ((pln->n - 1)/2);
177 pln->super.super.ops.add += 4 * ((pln->n - 1)/2);
178 pln->super.super.ops.other += 1; /* estimator hack for nop plans */
179
180 return &(pln->super.super);
181 }
182
183 /* constructor */
mksolver(void)184 static solver *mksolver(void)
185 {
186 static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 };
187 S *slv = MKSOLVER(S, &sadt);
188 return &(slv->super);
189 }
190
X(dft_r2hc_register)191 void X(dft_r2hc_register)(planner *p)
192 {
193 REGISTER_SOLVER(p, mksolver());
194 }
195