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 "rdft/rdft.h"
23 #include <stddef.h>
24
destroy(problem * ego_)25 static void destroy(problem *ego_)
26 {
27 problem_rdft *ego = (problem_rdft *) ego_;
28 #if !defined(STRUCT_HACK_C99) && !defined(STRUCT_HACK_KR)
29 X(ifree0)(ego->kind);
30 #endif
31 X(tensor_destroy2)(ego->vecsz, ego->sz);
32 X(ifree)(ego_);
33 }
34
kind_hash(md5 * m,const rdft_kind * kind,int rnk)35 static void kind_hash(md5 *m, const rdft_kind *kind, int rnk)
36 {
37 int i;
38 for (i = 0; i < rnk; ++i)
39 X(md5int)(m, kind[i]);
40 }
41
hash(const problem * p_,md5 * m)42 static void hash(const problem *p_, md5 *m)
43 {
44 const problem_rdft *p = (const problem_rdft *) p_;
45 X(md5puts)(m, "rdft");
46 X(md5int)(m, p->I == p->O);
47 kind_hash(m, p->kind, p->sz->rnk);
48 X(md5int)(m, X(ialignment_of)(p->I));
49 X(md5int)(m, X(ialignment_of)(p->O));
50 X(tensor_md5)(m, p->sz);
51 X(tensor_md5)(m, p->vecsz);
52 }
53
recur(const iodim * dims,int rnk,R * I)54 static void recur(const iodim *dims, int rnk, R *I)
55 {
56 if (rnk == RNK_MINFTY)
57 return;
58 else if (rnk == 0)
59 I[0] = K(0.0);
60 else if (rnk > 0) {
61 INT i, n = dims[0].n, is = dims[0].is;
62
63 if (rnk == 1) {
64 /* this case is redundant but faster */
65 for (i = 0; i < n; ++i)
66 I[i * is] = K(0.0);
67 } else {
68 for (i = 0; i < n; ++i)
69 recur(dims + 1, rnk - 1, I + i * is);
70 }
71 }
72 }
73
X(rdft_zerotens)74 void X(rdft_zerotens)(tensor *sz, R *I)
75 {
76 recur(sz->dims, sz->rnk, I);
77 }
78
79 #define KSTR_LEN 8
80
X(rdft_kind_str)81 const char *X(rdft_kind_str)(rdft_kind kind)
82 {
83 static const char kstr[][KSTR_LEN] = {
84 "r2hc", "r2hc01", "r2hc10", "r2hc11",
85 "hc2r", "hc2r01", "hc2r10", "hc2r11",
86 "dht",
87 "redft00", "redft01", "redft10", "redft11",
88 "rodft00", "rodft01", "rodft10", "rodft11"
89 };
90 A(kind >= 0 && kind < sizeof(kstr) / KSTR_LEN);
91 return kstr[kind];
92 }
93
print(const problem * ego_,printer * p)94 static void print(const problem *ego_, printer *p)
95 {
96 const problem_rdft *ego = (const problem_rdft *) ego_;
97 int i;
98 p->print(p, "(rdft %d %D %T %T",
99 X(ialignment_of)(ego->I),
100 (INT)(ego->O - ego->I),
101 ego->sz,
102 ego->vecsz);
103 for (i = 0; i < ego->sz->rnk; ++i)
104 p->print(p, " %d", (int)ego->kind[i]);
105 p->print(p, ")");
106 }
107
zero(const problem * ego_)108 static void zero(const problem *ego_)
109 {
110 const problem_rdft *ego = (const problem_rdft *) ego_;
111 tensor *sz = X(tensor_append)(ego->vecsz, ego->sz);
112 X(rdft_zerotens)(sz, UNTAINT(ego->I));
113 X(tensor_destroy)(sz);
114 }
115
116 static const problem_adt padt =
117 {
118 PROBLEM_RDFT,
119 hash,
120 zero,
121 print,
122 destroy
123 };
124
125 /* Dimensions of size 1 that are not REDFT/RODFT are no-ops and can be
126 eliminated. REDFT/RODFT unit dimensions often have factors of 2.0
127 and suchlike from normalization and phases, although in principle
128 these constant factors from different dimensions could be combined. */
nontrivial(const iodim * d,rdft_kind kind)129 static int nontrivial(const iodim *d, rdft_kind kind)
130 {
131 return (d->n > 1 || kind == R2HC11 || kind == HC2R11
132 || (REODFT_KINDP(kind) && kind != REDFT01 && kind != RODFT01));
133 }
134
X(mkproblem_rdft)135 problem *X(mkproblem_rdft)(const tensor *sz, const tensor *vecsz,
136 R *I, R *O, const rdft_kind *kind)
137 {
138 problem_rdft *ego;
139 int rnk = sz->rnk;
140 int i;
141
142 A(X(tensor_kosherp)(sz));
143 A(X(tensor_kosherp)(vecsz));
144 A(FINITE_RNK(sz->rnk));
145
146 if (UNTAINT(I) == UNTAINT(O))
147 I = O = JOIN_TAINT(I, O);
148
149 if (I == O && !X(tensor_inplace_locations)(sz, vecsz))
150 return X(mkproblem_unsolvable)();
151
152 for (i = rnk = 0; i < sz->rnk; ++i) {
153 A(sz->dims[i].n > 0);
154 if (nontrivial(sz->dims + i, kind[i]))
155 ++rnk;
156 }
157
158 #if defined(STRUCT_HACK_KR)
159 ego = (problem_rdft *) X(mkproblem)(sizeof(problem_rdft)
160 + sizeof(rdft_kind)
161 * (rnk > 0 ? rnk - 1u : 0u), &padt);
162 #elif defined(STRUCT_HACK_C99)
163 ego = (problem_rdft *) X(mkproblem)(sizeof(problem_rdft)
164 + sizeof(rdft_kind) * (unsigned)rnk, &padt);
165 #else
166 ego = (problem_rdft *) X(mkproblem)(sizeof(problem_rdft), &padt);
167 ego->kind = (rdft_kind *) MALLOC(sizeof(rdft_kind) * (unsigned)rnk, PROBLEMS);
168 #endif
169
170 /* do compression and sorting as in X(tensor_compress), but take
171 transform kind into account (sigh) */
172 ego->sz = X(mktensor)(rnk);
173 for (i = rnk = 0; i < sz->rnk; ++i) {
174 if (nontrivial(sz->dims + i, kind[i])) {
175 ego->kind[rnk] = kind[i];
176 ego->sz->dims[rnk++] = sz->dims[i];
177 }
178 }
179 for (i = 0; i + 1 < rnk; ++i) {
180 int j;
181 for (j = i + 1; j < rnk; ++j)
182 if (X(dimcmp)(ego->sz->dims + i, ego->sz->dims + j) > 0) {
183 iodim dswap;
184 rdft_kind kswap;
185 dswap = ego->sz->dims[i];
186 ego->sz->dims[i] = ego->sz->dims[j];
187 ego->sz->dims[j] = dswap;
188 kswap = ego->kind[i];
189 ego->kind[i] = ego->kind[j];
190 ego->kind[j] = kswap;
191 }
192 }
193
194 for (i = 0; i < rnk; ++i)
195 if (ego->sz->dims[i].n == 2 && (ego->kind[i] == REDFT00
196 || ego->kind[i] == DHT
197 || ego->kind[i] == HC2R))
198 ego->kind[i] = R2HC; /* size-2 transforms are equivalent */
199
200 ego->vecsz = X(tensor_compress_contiguous)(vecsz);
201 ego->I = I;
202 ego->O = O;
203
204 A(FINITE_RNK(ego->sz->rnk));
205
206 return &(ego->super);
207 }
208
209 /* Same as X(mkproblem_rdft), but also destroy input tensors. */
X(mkproblem_rdft_d)210 problem *X(mkproblem_rdft_d)(tensor *sz, tensor *vecsz,
211 R *I, R *O, const rdft_kind *kind)
212 {
213 problem *p = X(mkproblem_rdft)(sz, vecsz, I, O, kind);
214 X(tensor_destroy2)(vecsz, sz);
215 return p;
216 }
217
218 /* As above, but for rnk <= 1 only and takes a scalar kind parameter */
X(mkproblem_rdft_1)219 problem *X(mkproblem_rdft_1)(const tensor *sz, const tensor *vecsz,
220 R *I, R *O, rdft_kind kind)
221 {
222 A(sz->rnk <= 1);
223 return X(mkproblem_rdft)(sz, vecsz, I, O, &kind);
224 }
225
X(mkproblem_rdft_1_d)226 problem *X(mkproblem_rdft_1_d)(tensor *sz, tensor *vecsz,
227 R *I, R *O, rdft_kind kind)
228 {
229 A(sz->rnk <= 1);
230 return X(mkproblem_rdft_d)(sz, vecsz, I, O, &kind);
231 }
232
233 /* create a zero-dimensional problem */
X(mkproblem_rdft_0_d)234 problem *X(mkproblem_rdft_0_d)(tensor *vecsz, R *I, R *O)
235 {
236 return X(mkproblem_rdft_d)(X(mktensor_0d)(), vecsz, I, O,
237 (const rdft_kind *)0);
238 }
239