To check out this repository please hg clone the following URL, or open the URL using EasyMercurial or your preferred Mercurial client.

The primary repository for this project is hosted at https://github.com/sonic-visualiser/sv-dependency-builds .
This repository is a read-only copy which is updated automatically every hour.

Statistics Download as Zip
| Branch: | Tag: | Revision:

root / src / fftw-3.3.8 / dft / ct.c @ 167:bd3cc4d1df30

History | View | Annotate | Download (6.15 KB)

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

    
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

    
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

    
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

    
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

    
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

    
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

    
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

    
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

    
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

    
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
}