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 / rank-geq2.c @ 167:bd3cc4d1df30

History | View | Annotate | Download (5.53 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
/* plans for DFT of rank >= 2 (multidimensional) */
23

    
24
#include "dft/dft.h"
25

    
26
typedef struct {
27
     solver super;
28
     int spltrnk;
29
     const int *buddies;
30
     size_t nbuddies;
31
} S;
32

    
33
typedef struct {
34
     plan_dft super;
35

    
36
     plan *cld1, *cld2;
37
     const S *solver;
38
} P;
39

    
40
/* Compute multi-dimensional DFT by applying the two cld plans
41
   (lower-rnk DFTs). */
42
static void apply(const plan *ego_, R *ri, R *ii, R *ro, R *io)
43
{
44
     const P *ego = (const P *) ego_;
45
     plan_dft *cld1, *cld2;
46

    
47
     cld1 = (plan_dft *) ego->cld1;
48
     cld1->apply(ego->cld1, ri, ii, ro, io);
49

    
50
     cld2 = (plan_dft *) ego->cld2;
51
     cld2->apply(ego->cld2, ro, io, ro, io);
52
}
53

    
54

    
55
static void awake(plan *ego_, enum wakefulness wakefulness)
56
{
57
     P *ego = (P *) ego_;
58
     X(plan_awake)(ego->cld1, wakefulness);
59
     X(plan_awake)(ego->cld2, wakefulness);
60
}
61

    
62
static void destroy(plan *ego_)
63
{
64
     P *ego = (P *) ego_;
65
     X(plan_destroy_internal)(ego->cld2);
66
     X(plan_destroy_internal)(ego->cld1);
67
}
68

    
69
static void print(const plan *ego_, printer *p)
70
{
71
     const P *ego = (const P *) ego_;
72
     const S *s = ego->solver;
73
     p->print(p, "(dft-rank>=2/%d%(%p%)%(%p%))",
74
              s->spltrnk, ego->cld1, ego->cld2);
75
}
76

    
77
static int picksplit(const S *ego, const tensor *sz, int *rp)
78
{
79
     A(sz->rnk > 1); /* cannot split rnk <= 1 */
80
     if (!X(pickdim)(ego->spltrnk, ego->buddies, ego->nbuddies, sz, 1, rp))
81
          return 0;
82
     *rp += 1; /* convert from dim. index to rank */
83
     if (*rp >= sz->rnk) /* split must reduce rank */
84
          return 0;
85
     return 1;
86
}
87

    
88
static int applicable0(const solver *ego_, const problem *p_, int *rp)
89
{
90
     const problem_dft *p = (const problem_dft *) p_;
91
     const S *ego = (const S *)ego_;
92
     return (1
93
             && FINITE_RNK(p->sz->rnk) && FINITE_RNK(p->vecsz->rnk)
94
             && p->sz->rnk >= 2
95
             && picksplit(ego, p->sz, rp)
96
          );
97
}
98

    
99
/* TODO: revise this. */
100
static int applicable(const solver *ego_, const problem *p_, 
101
                      const planner *plnr, int *rp)
102
{
103
     const S *ego = (const S *)ego_;
104
     const problem_dft *p = (const problem_dft *) p_;
105

    
106
     if (!applicable0(ego_, p_, rp)) return 0;
107

    
108
     if (NO_RANK_SPLITSP(plnr) && (ego->spltrnk != ego->buddies[0])) return 0;
109

    
110
     /* Heuristic: if the vector stride is greater than the transform
111
        sz, don't use (prefer to do the vector loop first with a
112
        vrank-geq1 plan). */
113
     if (NO_UGLYP(plnr))
114
          if (p->vecsz->rnk > 0 &&
115
              X(tensor_min_stride)(p->vecsz) > X(tensor_max_index)(p->sz))
116
               return 0;
117

    
118
     return 1;
119
}
120

    
121
static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
122
{
123
     const S *ego = (const S *) ego_;
124
     const problem_dft *p;
125
     P *pln;
126
     plan *cld1 = 0, *cld2 = 0;
127
     tensor *sz1, *sz2, *vecszi, *sz2i;
128
     int spltrnk;
129

    
130
     static const plan_adt padt = {
131
          X(dft_solve), awake, print, destroy
132
     };
133

    
134
     if (!applicable(ego_, p_, plnr, &spltrnk))
135
          return (plan *) 0;
136

    
137
     p = (const problem_dft *) p_;
138
     X(tensor_split)(p->sz, &sz1, spltrnk, &sz2);
139
     vecszi = X(tensor_copy_inplace)(p->vecsz, INPLACE_OS);
140
     sz2i = X(tensor_copy_inplace)(sz2, INPLACE_OS);
141

    
142
     cld1 = X(mkplan_d)(plnr, 
143
                        X(mkproblem_dft_d)(X(tensor_copy)(sz2),
144
                                           X(tensor_append)(p->vecsz, sz1),
145
                                           p->ri, p->ii, p->ro, p->io));
146
     if (!cld1) goto nada;
147

    
148
     cld2 = X(mkplan_d)(plnr, 
149
                        X(mkproblem_dft_d)(
150
                             X(tensor_copy_inplace)(sz1, INPLACE_OS),
151
                             X(tensor_append)(vecszi, sz2i),
152
                             p->ro, p->io, p->ro, p->io));
153
     if (!cld2) goto nada;
154

    
155
     pln = MKPLAN_DFT(P, &padt, apply);
156

    
157
     pln->cld1 = cld1;
158
     pln->cld2 = cld2;
159

    
160
     pln->solver = ego;
161
     X(ops_add)(&cld1->ops, &cld2->ops, &pln->super.super.ops);
162

    
163
     X(tensor_destroy4)(sz1, sz2, vecszi, sz2i);
164

    
165
     return &(pln->super.super);
166

    
167
 nada:
168
     X(plan_destroy_internal)(cld2);
169
     X(plan_destroy_internal)(cld1);
170
     X(tensor_destroy4)(sz1, sz2, vecszi, sz2i);
171
     return (plan *) 0;
172
}
173

    
174
static solver *mksolver(int spltrnk, const int *buddies, size_t nbuddies)
175
{
176
     static const solver_adt sadt = { PROBLEM_DFT, mkplan, 0 };
177
     S *slv = MKSOLVER(S, &sadt);
178
     slv->spltrnk = spltrnk;
179
     slv->buddies = buddies;
180
     slv->nbuddies = nbuddies;
181
     return &(slv->super);
182
}
183

    
184
void X(dft_rank_geq2_register)(planner *p)
185
{
186
     static const int buddies[] = { 1, 0, -2 };
187
     size_t i;
188

    
189
     for (i = 0; i < NELEM(buddies); ++i)
190
          REGISTER_SOLVER(p, mksolver(buddies[i], buddies, NELEM(buddies)));
191

    
192
     /* FIXME:
193

194
        Should we try more buddies? 
195

196
        Another possible variant is to swap cld1 and cld2 (or rather,
197
        to swap their problems; they are not interchangeable because
198
        cld2 must be in-place).  In past versions of FFTW, however, I
199
        seem to recall that such rearrangements have made little or no
200
        difference.
201
     */
202
}